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by 

Anthony  V.  Fiacco 
Abolfazl  Ghaemi 

1.  Introduction 

This  is  a  manual  intended  to  facilitate  the  use  of  the 
computer  program  SENSUMT.  This  program  is  capable  of  solving  general 
parametric  nonlinear  programs  and  conducting  a  sensitivity  analysis 
using  the  Sequential  Unconstrained  Minimization  Technique  [10].  It  has 
recently  been  modified  to  include  the  calculation  of  piecewise  linear 
parametric  bounds  on  the  optimal  value  function  of  classes  of  nonlinear 
programming  problems  for  which  this  function  is  convex  or  concave.  In 
addition,  it  can  calculate  similar  bounds  on  a  class  of  separable  non- 
convex  right-hand  side  programs,  once  appropriate  convex  overestimating 
and  underestimating  programs  are  formulated. 

SENSUMT  is  the  outgrowth  of  a  routine  [3]  which  first  implemented 
a  penalty  function  sensitivity  approach.  This  was  motivated  by  results 
given  by  Fiacco  and  McCormick  [10]  for  a  particular  class  of  perturba¬ 
tions.  The  algorithmic  procedure  for  this  routine  was  suggested  by 
Fiacco.  It  was  later  refined  and  recoded  by  My lander  [12]  making  use  of 


several  subroutines  from  the  SUMT-Version  4  computer  program  coded  by 
My lander.  Holmes,  and  McCormick  [13].  These  routines  were  subsequently 
completely  integrated  with  SUMT-Version  4  by  Armacost  and  My lander  [7]. 
Armacost  [i]  further  revised  the  program  to  implement  the  generalized 
version  of  the  sensitivity  theory  developed  by  Fiacco  [9]  for  general 
parameter  perturbations.  This  routine,  again  following  a  procedure  sug¬ 
gested  by  Fiacco,  was  recently  modified  by  Ghaemi  to  conduct  piecewise 
linear  parametric  bound  calculation  on  the  optimal  value  function  of 
certain  classes  of  nonlinear  programming  problems  subject  to  large  para¬ 
metric  changes  [11]. 

This  latter  version  of  the  model,  designated  "SENSUMT,"  is  com¬ 
piled  at  the  Center  for  Academic  and  Administrative  Computing  of  The 
Ceorge  Washington  University. 

The  various  sections  of  this  manual  are  arranged  as  follows. 
Section  2  reviews  the  basic  sensitivity  results  and  presents  the  steps 
of  the  algorithm  implementing  sensitivity  calculation  via  a  sequential 
unconstrained  minimization  algorithm. 

Section  3  briefly  reviews  the  procedure  for  calculating  piece- 
wise  linear  parametric  bounds  on  the  optimal  value  function  of  certain 
classes  of  parametric  nonlinear  programming  problems.  It  presents  the 
relevant  algorithms  that  have  been  implemented. 

Section  4  gives  the  details  regarding  the  procedure  for  coding 
the  problems  for  solution,  sensitivity  analysis,  and  bound  calculation 
with  SENSUMT. 

Section  5  provides  the  codes  and  corresponding  computer  solution 
for  two  illustrative  examples  that  are  taken  from  I'll].  For  clarity,  an 
annotated  computer  listing  of  the  input  and  output  for  the  first  example 
is  also  provided. 

Section  6  gives  a  description  and  listing  of  the  various  sub¬ 
routines  comprising  SENSUMT.  This  section  is  included  to  make  the 


manual  self-contained.  With  the  exception  of  the  subroutines  BOUND, 
PERT,  and  TRANS,  and  the  new  version  of  program  MAIN,  the  material  in 
this  section  with  a  few  minor  modifications  has  been  taken  from  [1]  and 

U3]. 


2.  Sensitivity  Analysis  in  Nonlinear  Programming 
2.1  Basic  Sensitivity  Results 

The  parametric  mathematical  programming  problem  considered  by 
Fiacco  [9]  is  of  the  following  general  form: 

minimize  f(x,£) 
x£En 

P(£) 

subject  to  g^x.e)  >  0  ,  i=l,...,m  , 

hj(x,e)  =  0  ,  j=l,...,p  , 

where  x  is  the  usual  vector  of  variables  and  t  is  a  k-component  vec¬ 
tor  of  numbers  called  "parameters."  It  is  desired  to  analyze  the  behav¬ 
ior  of  a  solution  vector  x(e)  and  the  optimal  solution  value  f*(r)  = 
f[x(e),e]  near  some  given  value  of  e  .  Without  loss  of  generality, 
assume  that  the  parameter  vector  of  interest  is  e=0  . 

The  Lagrangian  for  Problem  P(e)  is  defined  as 

ra  P 

L(x,u,w,e)  =  f(x,£)  -  l  u  g  (x,e)  +  l  w  h  (x,e)  .  (2.1) 

1=1  11  j=l  3  2 

The  sensitivity  results  are  based  on  the  following  four  assump¬ 
tions  : 

A1  —  The  functions  defining  Problem  P(e)  are  twice  contin¬ 
uously  differentiable  in  (x,e)  in  a  neighborhood  of 
(x*,0)  . 

A2  —  The  second  order  sufficient  conditions  for  a  local 

minimum  of  Problem  P(0)  hold  at  x*  with  associated 
Lagrange  multipliers  y*  and  w*  . 
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A3  —  The  gradients  (x*>0)  >  for  all  i  such  that 

gi(x*,0)  =  0  ,  and  V^x^O)  ,  j-i,...,p  are 
linearly  independent. 

A4  —  Strict  complementary  slackness  holds  at  x*  when 
e=0  [i.e.,  u* > 0  for  all  i  such  that 

8i(x*,0)  =  0  ]. 

Under  the  above  assumptions,  Fiacco  [9]  established  the  following  gener¬ 
alization  of  Theorem  6  in  [101. 

Lemma  2.1  [local  characterization  of  a  Kuhn-Tucker  triple]:  If 
assumptions  Al,  A2,  A3,  and  A4  hold  for  Problem  F(e)  at  (x*,0)  ,  then 

(a)  x*  is  a  local  isolated  minimizing  point  of  Problem  P(0) 
and  the  associated  Lagrange  multipliers  u*  and  w* 
are  unique; 

(b)  for  e  in  a  neighborhood  of  0,  there  exists  a  unique 

once  continuously  differentiable  vector  function  y(e)  = 

T 

(x(e) ,u(e) ,w(e))  satisfying  the  second  order  sufficient 

conditions  for  a  local  minimum  of  Problem  P(e)  such  that 
T 

y(0)  =  (x*,u*,w*)  =  y*  ,  and  hence,  x(e)  is  a  locally 

unique,  local  minimum  of  Problem  P(e''  with  associated 
unique  Lagrance  multipliers  u(e)  and  w(e)  ;  and 

(c)  for  e  near  0,  the  set  of  binding  inequalities  is  un¬ 
changed,  strict  complementary  slackness  continues  to 
hold,  and  the  binding  constraint  gradients  are  linearly 
independent  at  x(e)  . 

(d)  (Armacost  and  Fiacco  [3]),  for  e  near  0,  the  gradient  . 
of  the  optimal  value  function  is 

V*(C)  =  VEL(y(E)’e)  ’  (2.2) 

(e)  which  also  means  that,  for  e  near  0,  the  Hessian  of  the 
optimal  value  function  is 

V^f*(e)  -  V*L(y(e),e)  .  (2.3) 

-  4  - 


4 


j 

e 
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The  above  results  provide  a  characterization  of  a  local  solution 
of  Problem  P(£)  and  its  associated  optimal  Lagrange  multipliers  near 
£=0  .  They  show  that  the  Kuhn-Tucker  triple  y(e)  is  unique  and  well 
behaved,  under  the  given  conditions.  Since  y(e)  is  once  differen¬ 
tiable,  the  partial  derivatives  of  the  components  of  y(e)  are  well  de¬ 
fined.  This  fact  and  assumption  A1  also  mean  that  the  functions  defin¬ 
ing  Problem  P(e)  are  once  continuously  differentiable  functions  of  e 
along  the  "solution  trajectory"  x(e)  near  £=0  ,  and  the  Lagrangian  is 
a  once  continuously  differentiable  function  of  £  along  the  "Kuhn- 
Tucker  point  trajectory."  The  above  results  constitute  the  structure 
for  numerous  developments  and  extensions,  many  of  which  have  been  estab¬ 
lished  by  Fiacco  {9]  and  Armacost  and  Fiacco  [2  -  6]. 


The  realization  of  this  theorem  for  the  parametric  right  hand 
side  problem  of  special  interest  in  the  present  study  is  treated  in  de¬ 
tail  by  Armacost  and  Fiacco  [4],  The  parametric  right  hand  side  problem 
is  the  following  important  realization  of  Pt£): 


minimize  f(x) 

subject  to  g^x)  >  ££  ,  i-1, . . .  ,m  , 


R(e) 


(x)  ^j+m  ’  J  1»*«**P  • 


The  Lagrangian  for  R(e)  is 


ux,u,»)  E  fW  -  “i(®iCx>-ei) +  • 

As  is  evident  from  the  Lagrangian,  the  results  (d)  and  (e)  of  Lemma  2.1 
for  Problem  R(£)  simplify  respectively  to 


and 


(d’)  Vef*(E) 


(e’)  V*f*(£) 


Veu(c) 

-Vew(e) 


(2.4) 


(2.5) 
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Fiacco  [9]  has  shown  that  the  class  of  algorithms  based  on  twice 
continuously  differentiable  penalty  functions  (specifically,  using  the 
logarithmic-quadratic  loss  penalty  function)  can  be  used  to  estimate 
y(e)  and  its  derivatives  in  a  neighborhood  of  c.*0  ,  for  the  general 
problem  P(e).  Minimization  of  the  penalty  function  with  penalty  param¬ 
eter  r  yields  a  solution  of  a  perturbation  of  the  Kuhn-Tucker  system 
in  a  neighborhood  of  (e,r)  =  (0,0)  .  Armacost  and  Fiacco  [3]  define  an 
optimal  value  penalty  function  and  obtain  first-  and  second-order  sensi¬ 
tivity  estimates  which  converge  to  the  corresponding  sensitivities  for 
the  optimal  value  function  for  Problem  P(e). 

The  logarithmic-quadratic  penalty  function  is 

m  P  „ 

W(x,e,r)  -  f(x,e)  -  r  l  £n  g  (x,e)  +  (J/2r)  l  h  (x,e)  .  (2.6) 

i=l  1  j*=l  J 

Lemma  2.2  (Fiacco  [9,  Theorem  3.1]):  If  the  assumptions  A1 
through  A4  hold,  then  in  a  neighborhood  of  (e,r)  =  (0,0)  there  exists 
a  unique  once  continuously  differentiable  vector  function  y(£,r)  = 
[x(e,r),u(e,r),w(e,r)]T  satisfying 

VxL(x,u,w,e)  =  0  , 

uigi(x,e)  =  r  ,  i=l,...,m  ,  (2.7) 

yx*c)  -  v  -  j-1 . ?  • 

with  y(0,0)  =  (x*,u*,w*)  ,  and  such  that  for  any  (e,r)  near  (0,0) 
and  r>0  ,  x(e,r)  is  a  locally  unique  unconstrained  local  minimizing 
point  of  W(x,£,r)  ,  with  g^[x(e,r) ,e]  >  0  ,  i*l,...,m  ,  and 

2 

VxW[x(c,r) ,e,r]  positive  definite. 

The  relevance  of  Equations  (2.7)  is  the  fact  that,  under  the 
given  conditions,  when  r“0  ,  they  are  necessary  conditions  that  B*,ust 
hold  at  a  local  solution  of  P(0)  and,  with  r>0  ,  they  are  necessary 
conditions  for  an  unconstrained  minimum  of  W(x,e,r)  .  The  latter  fact 
can  be  made  obvious  by  solving  for  u^  and  Wj  in  (2.7)  and  obtaining 
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VxL(x,u,w,e) 


V  f  -  yu  V  g,  +  vw  v  h. 

X  ^  i  xei  L  j  X  j 

Vxf  -  r  I(l/gi)  V.,gi  +  (1/r)  yhjVxhj 
VxW(x,e,r)  . 


Thus,  if  y(e,r)  is  a  solution  of  (2.7),  tnen 


VxW[x(£,r),e,r] 


VxL[x(e,r),u(e,r),w(e,r),  ]  =  0  . 


(2.8) 


(2.9) 


This  explicit  connection  between  the  optimality  conditions  of  lo¬ 
cal  solutions  of  P(e)  and  unconstrained  minima  of  W(x,£,r)  makes  it 
possible  to  approximate  information  characterizing  a  local  solution  of 
P(e)  by  algorithmic  calculations  associated  with  utilizing  W(x,£,r)  to 
solve  P(£).  In  particular,  differentiating  (2.9)  with  respect  to  £ 
yields 


V^W  V  x  +  V 

X  £ 


2 

£X 


W 


0  , 


and  using  the  fact  that 
Lemma  2.2)  yields 


2 

V  W  is  positive  definite  (a  conclusion  of 


2-1  2 

V  x  =  -v  W  1  W 

£  X  EX 


(2.10) 


evaluated,  of  course,  at  x(£,r)  .  Given  Vex(E,r)  ,  the  derivatives  of 

the  multipliers,  V  u. (£,r)  and  V  w, (£,r)  ,  can  then  be  calculated  by 
£  i  £  j 

differentiating  the  last  two  systems  of  equations  of  (2.7)  at  x(£,r) 
with  respect  to  £  . 


Lemma  2 ■ 3  (Fiacco  [9]):  For  £  in  a  neighborhood  of  e=0  ,  it 
follows  that 

(a)  lim  y(£,r)  =  y(£,0)  =  y(£)  ,  the  Kuhn-Tucker  triple 
T-*0+ 

characterized  in  conclusion  (b)  of  Lemma  2.1;  and 

(b)  lim+  V£y(£,r)  =  Vey(£,0)  -  V£y(£)  . 
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Under  the  conditions  of  Lemma  2.1,  x(e,r)  is  a  locally  unique 
minimizing  point  of  W(x,e,r)  .  Define  the  "optimal  value  penalty  func¬ 
tion"  as 

W*(e,r)  =  W(x(e,r) ,e,r)  .  (2.11) 

Armacost  and  Fiacco  [2,  Theorem  4  and  Corollary  4.1]  have  obtained  fur¬ 
ther  results  useful  for  estimating  the  first-  an<l  second-order  sensitiv¬ 
ity  of  the  optimal  value  function  f*(e)  of  Problem  P(e). 

Lemma  2.4  (Armacost  and  Fiacco  [3]):  If  the  assumptions  A1 
through  A4  hold  for  Problem  P(e),  then  in  a  neighborhood  of  e=0  , 

(a;  lim  W*(e,r)  =  f*(e)  ; 
r-*0 

(b)  V£W*(e,r)  =  V£L(y,e)  at  y  =  y(e,r)  ; 

(c)  lim.  V  W*(e,r)  =  V.f*(e)  ; 
r->0 

(d)  V^w*(e,r)  =  V^L(y(e,r),e)  ;  and 

(e)  lim  vV(c,r)  =  V*f*(c)  ; 

r-K)  L 

where  convergence  is  component  by  component  in  all  cases. 

Lemmas  2.1  through  2.4  enable  us  to  calculate  an  estimate  of 

2 

y(e)  ,  Vcy(e)  ,  V£f *(c)  ,  and  V  f*(e)  when  e  is  near  0  and  r  is 
near  0,  once  y(e,r)  is  available. 

In  the  next  section  we  briefly  present  the  algorithmic  implemen¬ 
tation  of  some  of  the  above  results, 

2.2  The  Algorithm  Implementing  the 
Sensitivity  Results 

The  penalty  function  algorithm  SUMT  estimates  the  solution  of  the 
general  mathematical  problem  P (e)  by  estimating  the  unconstrained  minima 
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of  the  penalty  function  W(x,£,r)  at  successively  decreasing  values  of 
the  penalty  function  parameter  r  >  0  .  Under  conditions  weaker  than 
those  assumed  here,  Fiacco  and  McCormick  [10]  have  shown  that  as  r  ap¬ 
proaches  zero,  the  sequence  of  the  unconstrained  minima  of  W(x,£,r) 
will  approach  a  solution  of  P(£).  Each  unconstrained  penalty  function 
minimization  is  thus  a  "subproblem"  associated  with  a  particular  value 
of  the  penalty  function  parameter  r  . 

The  successive  steps  of  the  algorithm  for  first  order  sensitivity 
analysis  of  the  Kuhn-Tucker  triple  with  respect  to  the  j th  parameter 
that  are  implemented  in  SENSUMT  [1]  are  listed  below.  Here  the  notation 
x  or  x(e)  denotes  the  estimate  of  a  solution  point  of  P(E)  calculated 
by  SUMT  for  a  given  value  of  the  penalty  function  parameter  r  ,  where 
e  denotes  the  value  of  the  problem  parameter  at  which  this  sensitivity 
is  estimated. 


Algorithm  2.1: 

2  “1  2  —  —  -*1 

Step  1.  Compute  a  representation  of  V  W  »  VSl(x,£,r)  by  L-U 

xx  2 

decomposition  using  the  SUMT  subroutines.  If  V^W  is 
not  positive  definite,  terminate  the  sensitivity  analysis. 

T 

Step  2.  Estimate  3(V  W  )/3 £.  using  the  central  differencing 

x  3 

formula 

B(VxWT)/3e.  -  (i)(VxW(x,e+Ae^,r)T  -  V.;W(x,i-Aej  ,r)T )  , 

where  A  is  the  differencing  interval  and  e^  is  the  jth 
unit  vector. 


Step  3.  Calculate 


3x(c)/3e .  =  -vV1  3(V  WT)/3c. 

J  x  x  J 


Step  4.  Estimate  V^g^(x,z)  and  V^h^Cx.e)  using 

3g1(x,e)/3eJ  -  (-^(g^x.i+Aej)  -  g^i^i-Ae^))  ,  and 
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Step  5 


Step  6 


Step  7. 


9h1(x,S)/9e.j  •»  (-^)(hi(x,e+Ae^)  -  h^x.e-Aej))  . 

Estimate  the  components  of  V£u (E)  for  1*1,..., m  using 
8u4 (i)/9£j  «  -(r/g1(x,i)2 j(Vxgi(x,S)9x(i)/9ej  +  9gi(x,e)/9e^  . 

Estimate  the  components  of  V£w (E)  for  i=l . .  using 

9w1(£)/9£j  w  (1/r)  (v^h^Cx.i) 3x(£)/9Ej  V  9hi(x,e) '9^)  . 

There  are  tx^o  methods  for  estimating  7£f*(£)  ,  the  first  using 

7f**VfVx+Vf  ,  with  V  x  obtained  from  Step  3,  and  the 

e  x  c  e  £ 

second  method  using  the  gradient  of  the  penalty  function  W  , 
or  equivalently,  the  Lagrangian  taken  with  respect  to  the  pa¬ 
rameters  [Lemma  2.1,  conclusion  (d);  Lemma  2.4,  conclusion 
(b)j.  Both  are  incorporated  in  the  computer  program  but  used 
for  different  purposes.  The  former  method  gives  the  most  ac¬ 
curate  estimate  of  V£f*(i)  and  is  summarized  in  Steps  7  and  8. 

Estimate  the  components  of  V£f(x,e)  using  the  central  dif¬ 
ferencing  formula 

9f(x,e)/9€j  «  (-^)(f (x.i+Ae^)  -  f(x,i-Aej))  . 

Calculate  an  estimate  of  the  components  of  V£f*(e)  using  the 
results  of  Steps  3  and  7  as 

9f*(e)/9e.  «  V  f (x,e) 9x(e)/9e  +  9f(x,e)/9e  . 

J  x  J  J 

The  second  method,  using  the  gradient  of  the  Lagrangian  to  es¬ 
timate  Vcf*(e)  ,  is  computationally  less  expensive  and  is  used 
to  obtain  rough  estimates  that  single  out  the  more  crucial  pa¬ 
rameters  for  further  analysis.  This  approximation  is  calcu¬ 
lated  as  follows. 
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Step  9.  Estimate  the  components  of  using  the  results  of  Steps 

4  and  7  as 


m 

3f*(e)/3e.  as  3f(x,e)/3e.  -  l  u.(e,r)3g  (x,e)/3e 
J  3  i=l  1 


p 

+  l  wi(c,r)dh^(x,c) /dc .  . 
i=l  3 


3.  Parametric  Optimal  Value  Bounds 

In  this  section  we  briefly  review  the  procedure  for  calculating 
piecewise  linear  parametric  upper  and  lower  bounds  on  the  optimal  value 
function  of  convex  problems  with  right  hand  side  perturbations  of  the 
form  CR(e)  .  It  is  important  to  note  that  the  given  technique  is  not 
only  applicable  to  problems  of  this  form,  but  also  to  any  class  of  prob¬ 
lems  which  have  convex  or  concave  optimal  value  functions.  In  addition, 
we  show  how  the  above  technique  can  be  exterdea  to  calculate  similar 
bounds  on  the  optimal  value  function  of  separable  nonconvex  right  hand 
side  perturbation  problems  SR(e)  ,  once  their  (computable)  over-  and 
underestimating  problems  are  available.  For  a  more  detailed  treatment, 
readers  should  refer  to  [11]. 


3.1  Parametric  Bounds  on  the  Optimal  Value  Function 
of  Convex  Right  Hand  Side  Problems  CR(e) 


Consider  the  following  right  hand  side  parametric  programming 
problem  R(e),  discussed  in  Section  2: 


minimize  f(x) 

subject  to  g^x)  >_  ,  i=l,m 

hj  (x)  =  e1  ,  j=ntfl,p 


R(e) 


where  xeE  ,  f  ,  g  ,  and  h  :  E  -*■  E  ,  and  C  and  eeEp  .  If  f(x) 
-g^x)  ,  i“l,m  ,  are  convex  and  (x)  ,  j=m+l,p  are  linear,  then 


and 
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the  problem  R(£)  is  convex  and  will  be  designated  by  CR(£).  It  is  well 
known  that  f*(e)  ,  the  optimal  value  function  of  the  problem  CR(e),  is 
a  convex  function  of  £  . 

The  convexity  of  the  optimal  value  function  f*(£)  of  the  prob¬ 
lem  CR(e)  enables  one  to  calculate  parametric  upper  and  lower  bounds  on 
this  function  when  any  of  the  problem  parameters  is  radically  perturbed. 
This  of  course  requires  the  solution  and  corresponding  optimal  value 
function  sensitivity  information  for  both  perturbed  and  unperturbed 
problems.  Before  delving  into  the  calculation  procedure  of  these 
bounds,  we  must  remind  ourselves  of  the  following  two  basic  properties 
of  convex  functions. 

(i)  Any  line  connecting  two  points  on  the  graph  of  a  convex 
function  does  not  underestimate  that  function  between 
the  points. 

(ii)  Any  line  tangent  to  the  graph  of  a  convex  function  does 
not  overestimate  that  function. 

These  two  properties  lend  themselves  in  a  natural  way  to  the  calculation 
of  parametric  bounds  on  the  optimal  value  function  of  the  problem  CR(e) 
under  large  perturbations  of  any  of  the  problem  parameters,  say  £^  . 


3.2  Algorithm  for  Calculation  of  Bounds 
on  f*(£)  of  the  Problem  CR(e) 

In  the  following  we  list  the  successive  steps  required  in  calcu¬ 
lation  of  bounds  on  f*(e)  ,  the  optimal  value  function  of  CR(e),  as  a 
function  of  ,  when  is  perturbed  from  to  E^ 

while  the  remaining  parameters  are  fixed  at  their  base  values. 

Algorithm  3.1: 

Step  1.  Solve  the  unperturbed  problem  and  obtain  f*(e.)  and 
df*(ei)/dt-i  at  e^  =  .  Note  that 
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df*(£i)  | 

i  Ui(Li)  • 

j— l,...,m 

(2.4) 

d£i  ~i 

|-wt(Ei)  ’ 

i=nri-l,p  . 

(2.5) 

Re-solve  the  problem  and  obtain 

£i  =  Ei2  ‘ 

f * (e^ )  and 

d£*(.c^)/dz  at 

Step  Z.  Derive  the  equation  of  the  line  passing  through  the  points 

(e.  =  Ea  »  f*(ei}  =  f*(Eii})  and  (ci  =  ln  ’  f*(Ei}  = 

f*(e12))  .  This  line  provides  a  parametric  upper  bound  f*^^) 

for  f*(e.)  as  a  function  of  £.  £  [c.,  ,  £.„]  . 
i  1  il  i2 


Step  Derive  the  equation  of  the  tangent  lines  to  f*(£^)  at  the 
above  two  points  with  the  slopes 


df*(£i) 


£i=£il 


and 


df*(£.) 
v  1 

d£ 


£.=£ 

l 


12 


calculated  in  Steps  1  and  2,  respectively.  These  two  lines 
provide  a  parametric  lower  bound  f^E^)  for  f*(ei)  as  a 
function  of  e  [e^  ,  e^l  . 


The  lines  obtained  in  Steps  3  and  4  form  a  triangle  which  en¬ 
closes  the  optimal  value  function  f*(E^)  over  the  given  range  of  . 

An  estimate  of  f*(£^)  can  also  be  made  by  fitting  a  convex  function 

that  passes  through  the  points  given  in  Step  3  and  having  the  corre¬ 
sponding  slopes  at  these  points  obtained  in  Steps  1  and  2. 

It  is  clear  that  the  fundamental  requirement  in  using  the  above 
algorithm  is  that  the  optimal  value  function  be  convex  and  differenti¬ 
able  at  least  at  and  =  e^  •  as  it  was  given  in  Section 

2,  Fiacco,  under  the  assumptions  A1  -  A4  of  Lemma  2.1,  has  established 
the  differentiability  of  the  optimal  value  function  in  a  neighborhood  of 
£-0  for  the  general  parametric  programming  problem  P(£).  Thus,  not 
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only  Problem  CR(e),  but  all  classes  of  parametric  nonlinear  programming 
problems  possessing  a  convex  optimal  value  function  and,  in  particular, 
meeting  assumptions  A1  -  A4  at  and  E  =  e^,  can  analyzed 

for  bounds  according  to  Algorithm  3.1. 

3.3  Parametric  Bounds  on  the  Optimal  Value 
Function  of  the  Problem  SR(e) 

In  this  section  we  show  how  Algorithm  3.1  can  be  extended  to  cal¬ 
culate  piecewise  linear  parametric  upper  and  lower  bounds  on  the  optimal 
value  function  of  a  separable  nonconvex  right  hand  side  perturbation 
problem,  given  the  over-  and  underestimating  problems. 


3.3.1  Separable  right  hand  side  perturbation 
problems . 

The  separable  nonconvex  problem  addressed  here  has  the  following 
structure: 

n 

minimize  f(x)  =  J  f  (x.) 


subjec 


:ct  to  g  (x)  =  l  g.,(x .)  >€  ,  i=l, . . . , 

1  j=l  J  3 


SR(E) 


Lj  ~  Xj  =  Uj  ’  ’ 

where  each  f.(x.)  and  g.,(x,)  is  a  function  of  a  single  variable  x. 

J  J  ij  j  j 

and  is  differentiable;  and  are  lower  and  upper  bounds  on  the 

variable  x^  ,  respectively.  Note  that  the  discussion  that  follows  is 
also  valid  when  Problem  SR(e)  involves  linear  equality  constraints.  Let 


u 

G  =  x  :  g  (x)  =  l  g.,(x  )  >  e.  ,  1=1,..., m 

j-1  J  J 

Bj  E  ^Xj  1  Lj  =  Xj  =  M  ’  ' 


X  B  , 

j-1  J 
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and 


S  =  G  O  B  . 


The  set  G  ,  and  therefore  set  S  ,  depends  on  £  .  Thus  the  problem 
SR(e)  will  read  as  follows: 


minimize  f(x)  =  £  f  (x.) 


j=l 


r t 


subject  to  x  e  B^G  =  S  C  En  . 


SR(e) 


3.3.2  Convex  underestimating  problem. 


Let  f.(x.)  be  the  convex  envelope  of  f.(x.) 

r  i  j 

gjj(Xj)  be  the  concave  envelope  of  g„(x^)  ,  j=l,.. 

over  the  interval  B.  .  Thus,  for  x.  £  B.  , 

J  11 


and 


(xj )  >  8ij(xj)  *  i=1 . . 


,  j=l , . . . ,m  and 
,n  ,  1 , . . .  j m 


,n  . 


Definition:  Define  the  problem  CSR(£)  as  follows: 


where 


minimize  f(x)  =  £  f.(x.) 

x£En  ~  j=l  ~3  3 

subject  to  xeGOb  5  S, 

n  . 

*  :  gi(x)  =  l  g  (x)  >  Ei  ,  i=l, . . . ,m  >  . 


1=1 


It  follows  easily  that  problem  CSR(e)  is  a  separable  convex  RHS  pertur¬ 
bation  problem  whose  optimal  value  function  underestimates  the  global 
optimal  value  function  of  the  problem  SR(£) . 

Thus  the  basic  ingredients  required  to  formulate  the  convex  un¬ 
derestimating  problem  CSR(e)  are  the  convex  envelope  of  each  term  in  the 
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objective  function  and  the  concave  envelope  of  each  term  in  the  con¬ 
straints  of  the  problem  SR(e)  over  the  rectangular  polytope  B  . 

In  the  next  section  we  will  discuss  the  formulation  of  a  convex 
overestimating  problem  CSR(e)  of  the  problem  Sk(e) . 


3.3.3  Convex  overestimating  problem. 


Let  fj(Xj)  be  a  convex  function  which  does  not  underestimate 
fj(xj)  .  j=l,...,n  »  and  be  a  concave  function  which  does  not 

overestimate  g^Cx^)  ,  j=l,...,n  ,  i=l,...,m  ,  over  the  interval  . 
Thus,  for  Xj  e  Bj  , 


and 


Cj'V  i  VV 

g^ j (Xj )  <  g^j (xj )  »  i— l,...,mj  j-1, . . . ,n 


Henceforth  in  discussion,  f.(x.)  will  be  abbreviated  as  convex 

3 

"nuf"  (nonunderestimating  function)  of  the  single  variable  function 

f.(x,)  ,  and  g..(x.)  will  be  abbreviated  as  concave  "nof"  (nonoveres- 
J  j  iJ  j 

timating  function)  of  g^Cx^)  over  the  interval  ,  i=l,...,m  , 

j— l,...,n. 


Definition:  Let  us  define  the  problem  CSR(e)  as  follows: 


where 


minimize  f(x)  =  \  f.(x  ) 

x  En  *  j=l  J 

subject  to  x  e  (G  D  B)  =  I  , 

n 

:  L(x)  =  I  !,,(*,)  i  ei  » 
j-1  3  2 


It  is  easy  to  show  that  problem  CSR(e)  is  a  separable  convex  RHS  problem 
whose  optimal  value  function  underestimates  the  global  optimal  value 
function  of  the  problem  SR(e). 
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Thus,  in  order  to  formulate  an  overestimating  convex  problem 

CSR(e)  of  the  problem  SR(e),  one  must  calculate  convex  nuf,  f.(x.)  , 

*J  J 

for  each  component  of  f(x)  and  concave  nof,  g^^(Xj)  ,  for  each  compo¬ 


nent  of  g^(x)  ,  i=l,...,m  ,  over  the  interval  .  It  must  be  noted 

that  convex  nuf  and  concave  nof  of  a  given  function,  unlike  its  enve¬ 
lopes,  are  not  unique.  Thus  the  overestimating  problem  C§R(e)  of  prob¬ 
lem  SR(e)  is  not  unique.  For  further  treatment  of  convex  nuf  and  con¬ 
cave  nof  of  single  variable  functions,  readers  may  refer  to  [11]. 


An  important  point  to  be  noted  here  is  that  formulation  of  the 
convex  problems  CSR(e)  and  CSR(e)  enable  one  to  calculate  parametric 
upper  and  lower  bounds  on  the  (global)  optimal  value  function  of  the 
separable  nonconvex  programming  problem  SR(e). 


3.4  Algorithm  for  Calculation  of  Bounds  on 
the  (Clobal)  Optimal  Value  Function  of 
the  Problem  SR(e) 

In  the  following  an  algorithm  is  presented  for  calculation  of  up¬ 
per  and  lower  bounds  on  the  (global)  optimal  value  function  of  the  prob¬ 
lem  SR(e)  as  a  function  of  each  RHS  parameter  ei  ,  where  e  ,  e^l 

Implicit  in  this  algorithm  is  the  assumption  that  the  overestimating 
problem  CSR(E)  has  a  nonempty  interior,  which  is  a  basic  requirement  of 
almost  any  nonlinear  programming  algorithm. 


Algorithm  3.2: 

Step  1:  Calculate  f^(Xj)  ,  a  convex  nuf  for  each  f^(x^)  ,  j=l,...,n  , 
and  g^^)  ,  a  concave  nof  for  each  g^j(x^)  ,  i=l,...,m  , 
j=l,...,n  ,  in  closed  interval  B^  ,  and  construct  the  corres¬ 
ponding  overestimating  problem  CSR(e) . 


Step  2:  Solve  the  problem  C§R(e)  at  and  *  ei2  and  ob¬ 


tain  f*(e^)  at  these  values  of 


-  17  - 


■ 


T-434 


Step  3:  Derive  the  equation  for  the  line  passing  through  the  points 

(ei =  » **<ei)  =  and  (ei =  ^12  *  i*(ei}  ■ 

f*(E  2))  •  This  line  provides  a  parametric  upper  bound  f^e^) 

on  the  (global)  optimal  value  function  f*(e^)  of  the  problem 

SR(£)  as  a  function  of  over  the  interval  [e^  ,  ♦ 


Step  4:  Calculate  f^(Xj)  »  t*ie  convex  envelope  for  each  f^(Xj)  , 
i=l,...,n  and  g_  (xj )  »  t^ie  concave  envelope  for  each 
g^^Xj)  ,  i=l,...,m  ,  j=l,...,n  in  closed  interval  ,  and 
construct  the  corresponding  underestimating  problem  CSR(e) . 


Step  5: 


Solve  the  problem  CSR(e)  at  £^  =  £^  and  £^  =  and  ob¬ 

tain  f*(ei>  and  u^(e^)  =  df*(e^)/dei  at  these  values  of  e 


i  * 


Step  6:  Derive  the  equation  for  the  tangent  lines  to  f*(£^)  at  = 
c  and  *  e  ,  with  respective  slopes  of  ui(^il)  and 
ui(?-  2>  derived  in  Step  5.  These  two  lines  over  the  interval 
[e±l  ,  ei2]  jointly  provide  a  parametric  lower  bound  f*(e^) 
on  the  (global)  optimal  value  function  f*(e^)  of  the  problem 
SR(e)  as  a  function  of  , 


The  implicit  assumption  in  Steps  5  and  6  of  this  algorithm  is 
that  assumptions  A1  -  A4  of  Lemma  2.1  hold  at  the  solution  points  when 
£ ^  takes  values  of  and  • 

Before  delving  into  the  next  section,  we  must  point  out  that  if 
over-  and  underestimating  problems  with  concave  optimal  value  functions 
could  be  formulated,  then  the  algorithm  above,  with  minor  alterations, 
can  be  used  to  calculate  the  desired  bounds.  In  this  instance  Step  6, 
when  applied  to  the  overestimating  problem,  will  provide  an  upper  bound 


-  18  - 


T-434 


while  Step  3,  when  applied  to  the  underestimating  problem,  will  provide 
a  lower  bound  on  the  optimal  value  function.  For  calculation  of  enve¬ 
lopes,  convex  nuf  and  concave  nof  of  single  variable  functions,  readers 
may  refer  to  [11].  The  entries  in  Table  1  summarize  the  concept  of  con¬ 
vex  nuf  and  concave  nof  for  a  variety  of  single  variable  functions. 

4.  SENSUMT  Input  Specifications 

The  coding  procedure  for  nonlinear  programming  problems  for  solu¬ 
tion,  sensitivity  analysis,  and  bound  calculation  using  SENSUMT  closely 
follows  the  coding  procedure  for  SUMT-Version  4  [13],  In  order  to  code 
a  given  problem  for  the  above  calculations,  one  must  have  (i)  a  user- 
supplied  subroutine,  and  (ii)  user's  information  cards. 

4.1  User-supplied  Subroutines 

User-supplied  subroutines  will  generally  consist  of  four  sub¬ 
routines  called  READIN,  RESTNT,  GRAD1  and  MATRIX.  These  subroutines  are 
the  only  subroutines  for  SENSUMT  that  are  problem  dependent.  The  com¬ 
munication  between  these  subroutines  and  the  rest  of  the  subroutines  of 
SENSUMT  is  mostly  through  the  COMMON  blocks,  but  partly  via  their  argu¬ 
ments.  Each  of  the  user-supplied  subroutines  must  contain  the  double 
precision  card 

IMPLICIT  REAL*8(A-H,0-Z) . 

The  subroutines  RESTNT,  GRAD1,  and  MATRIX  must  contain  the  COMMON  card 

COMMON/ SHARE/X (45),  DEL(45),  A(45,45),  N,  M,  MN,  NP1,  NM1 . 

Also,  depending  on  the  problem  being  solved,  the  following  COMMON  card 
may  be  required  in  some  or  all  of  the  user  subroutines: 

COMMON/ SEN/PAR (45) ,  DPAR(45),  NPAR,  ISENS. 

When  bound  analysis  is  desired,  subroutine  READIN  must  contain  the  COMMON 
card 

COMMON/ABG/LY ,  LZ,  PER(45). 
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In  the  following  section  we  discuss  the  purpose  and  coding  procedure  of 
each  of  these  user-supplied  subroutines.  The  general  description  and 
listing  of  the  rest  of  the  subroutines  comprising  SENSUMT  will  be  given 
in  Seccion  6. 

4.1.1  READ IN 

This  subroutine  allows  the  user  to  read  in  the  data  which  is 
necessary  for  evaluation  of  the  problem  functions  and  convey  it,  through 
the  COMMON  blocks,  to  the  related  subroutines.  The  problem  parameters 
for  which  sensitivity  and  bound  calculation  are  desired  must  also  be  de¬ 
fined  and  initialized  in  this  subroutine.  The  initial  values  of  the 
parameters  must  be  placed  in  array  PAR.  The  corresponding  perturbation 
values  for  bound  calculation  must  be  placed  in  array  PER.  NPAR  in  this 
subroutine  must  be  initialized  to  the  total  number  of  problem  parameters 
on  which  sensitivity  is  desired.  Notice  that  the  presence  of  array  PER 
in  this  subroutine  will  trigger  SENSUMT  to  calculate  the  desired  bound 
calculation.  READIN  may  also  be  used  to  read  in  and  print  out  any  per¬ 
tinent  information  of  the  user's  choice,  such  as  problem  title,  date, 
and  so  on.  This  subroutine  is  called  only  once  for  each  problem  being 
solved.  For  clarity,  see  Example  1  and  its  code  in  Figure  2.  Any  input 
data  read  by  subroutine  READIN  must  be  placed  immediately  after  the 
first  option  card. 

4.1.2  RESTNT  (IN.VAL) 

This  subroutine  is  used  to  read  in  the  problem  functions.  For 
IN=0,  set  VAL=f(x)  and  for  IN^O  set 


VAL  =  8IN(x)  ’ 

IN=1,2, . . . .m 

VAL  =  hJN(x)  , 

IN=m+l, . . . ,m+p  . 

For  clarity  see  Example  1  and  Figure  2. 
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4.1.3  GRAD1  (IN) 

This  subroutine  is  used  to  read  in  the  gradients  of  the  problem 
functions.  For  IN=0  set 

DEL(J)  =  f(x)/3xJ  ,  J=l,2,...,n 

and  for  IN^O  set 

!9gIN(x)/3xJ  ,  J=l,2,...,n 

. 

3hIN(x)/3JcJ  ,  J=l,2,...,n 

For  an  illustration,  see  Example  1  and  Figure  2. 

SENSUMT  can  internally  compute  numerical  approximations  for  some 
or  all  of  the  gradients.  To  utilize  this  option,  instead  of  coding  the 
gradients  code  the  statements 

CALL  DIFFl(IN) 

RETURN 

for  any  desired  value  of  the  variable  IN  (the  index  of  the  problem  func¬ 
tions,  with  IN=0  corresponding  to  the  objective  function  f(x)  ).  Be¬ 
cause  of  the  computational  effort  involved  in  numerical  differencing, 
this  option  may  be  prohibitive  for  large  problems. 

4.1.4  MATRIX  (IN,K) 

This  subroutine  reads  in  the  upper  triangle  and  diagonal  elements 
of  the  Hessian  matrix  of  the  problem  functions.  For  IN=0,  set 

A(I,J)  =  32f (x)/3xI3Xj  ,  1=1,2,.. .,n;  J=l,2,...,n;  and  l£J  . 

For  IN^O  set 

32glN(x)/3Xi3xj  ,  IN-1,2,. ..,m 

A(I,J)  -  j  » 

3  hIN (x) / 3Xj  3Xj  ,  IN=m+l , . . . , p 

where  1=1,2,. ..,n;  J=l,2,...,n;  and  I  <  J  .  Before  the  call  is  made  to 
this  subroutine,  all  entries  of  the  matrix  A(I,J)  are  set  to  zero. 
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The  second  argument  K  of  the  subroutine  MATRIX  is  provided  so 
that  the  user  may  communicate  to  SENSUMT  that  the  Hessian  of  a  problem 
function  is  identically  zero.  Set  K=1  if  the  second  partial  deriva¬ 
tives  of  the  IN th  constraint  are  zero  (IN=0  corresponds  to  the  objective 
function).  For  an  illustration,  see  Example  1  and  Figure  2. 

SENSUMT  can  internally  compute  numerical  approximations  for  some 
or  all  of  the  Hessian  matrices.  To  utilize  this  option,  instead  of 
coding  the  Hessian  matrices,  code  the  statements 

CALL  DIFF2 (IN) 

RETURN 

for  any  desired  value  of  the  variable  IN  (the  index  of  the  problem  func¬ 
tion,  with  IN=0  designating  the  objective  function).  Although  this  op¬ 
tion  can  he  used  in  calculating  an  optimal  solution ,  the  current  sensi¬ 
tivity  routines  require  explicit  coding  of  the  closed  form  of  the  Hes¬ 
sian  matrices.  Thus,  this  option  cannot  be  used  for  sensitivity  and 
bound  calculations.  Again,  because  of  the  computational  effort  involved 
in  numerical  differencing,  use  of  this  option  is  not  advised  for  large 
problems. 

4.2  User's  Information  Cards 

The  other  inputs  required  by  SENSUMT  are  the  user's  information 
cards ,  i . e . , 

•  PARAMETER  CARD 

•  INITIAL  VECTOR  CARD(S) 

•  FIRST  OPTION  CARD 

•  TOLERANCE  CARD 

•  SECOND  OPTION  CARD. 

In  the  following  we  specify  the  type  of  information,  along  with  the  cor¬ 
responding  format,  name,  and  description,  which  must  be  provided  by  the 
user  to  SENSUMT  via  the  above  cards.  All  of  this  information  is  read  by 
the  program  MAIN,  which  sets  up  the  SENSUMT  algorithm  to  solve  and 
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analyze  the  problem  under  study.  Tables  2,  3,  and  4,  which  follow,  are 
taken  from  [13]. 


4.2.1  Parameter  card. 


Input  specifications  of  the  parameter  card  are  shown  in  Table  2. 


TABLE  2 

PARAMETER  CARD 


Columns 

Format 

Name 

Use 

01-12 

E12.0 

epsi  (e) 

Tolerance  used  to  decide  if  an  un¬ 
constrained  minimum  has  been 
achieved  for  each  subproblem  [see 
Option  9] 

13-24 

E12.0 

RH0IN  (rx) 

Possible  initial  value  of  r  (often 
set  at  1.0)  [see  Option  1] 

25-36 

E12.0 

THETAO  (0Q) 

Tolerance  used  tc  decide  if  the 
solution  to  the  NLP  problem  P(e) 
has  been  approximated  [see  Option 

5] 

37-48 

E12.0 

RAT 10  (c) 

Parameter  (>1)  used  to  compute 
consecutive  values  of  r;  r  .  = 
r^/c  (often  set  at  16.0)  1 

49-60 

E12.0 

TMMAX 

Maximum  amount  of  time  for  solving 
problem  (in  seconds) 

61-64 

14 

M 

Number  (integer)  of  nontrivial 
constraints  [see  Option  2] 

(M+-MZ)  <  200* 

65-68 

14 

N 

Number  (integer)  of  variables,  N  < 

45 

69-72 

14 

MZ 

. 

Number  (integer)  of  equality  con¬ 
straints 

*'L'he  limits  on  M+MZ  and  N  are  governed  by  the  size  of  the  arrays  in 
SKNSUMT. 

4.2.2  Initial  vector  card(s). 

The  cards  designating  the  initial  starting  point  immediately  fol¬ 
low  the  parameter  card.  There  are  six  components  per  card,  requiring 
n/6  cards  for  the  initial  vector.  Each  card  has  the  format  6E12.6  . 


-  24  - 


4.2.3  First  option  card. 


The  input  specifications  for  the  first  option  card  are  shown  in 
Table  3.  This  card  must  immediately  follow  the  user-supplied  subrou¬ 
tines. 

4.2.4  Tolerance  card. 

The  input  specifications  for  the  tolerance  card  are  given  in 
Table  4.  This  card  must  immediately  follow  the  first  option  card. 

4.2.5  Second  option  card. 

The  input  specifications  for  the  second  option  card,  which  im¬ 
mediately  follows  the  tolerance  card,  are  given  in  Table  5. 

Figure  1  shows  the  arrangement  of  the  data  together  with  JOB  and 
JCL  cards  in  a  coded  deck. 
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Option 


Column 


Meaning 


(normally 
set  to  3) 


The  value  for  r  is  made  by  finding 
an  approximate  solution 
min{ [Vw(x° ,r) [ ?2W(x° ,r)  ]"  !Vw(x° ,r) ] 
which  is  a  good  approximation  only 
when  x°  is  close  to  the  boundary 
(i.e.,  for  some  i,  g^(x)  is  close  to 
zero)  or  when  V2f(x°)=0  and  when 
MZ=0 . 

Tj  is  given  by  formula  8.65  [10,  p. 
191]  (can  only  be  used  when  MZ=0) . 
rj  =  RH0IN  (see  parameter  card). 

The  requirements  (trivial  con¬ 
straints)  that  x-l  J>  0  for  i=l,...,n 
are  to  be  automatically  included  in 
the  problem. 

The  only  constraints  on  the  problem 
are  those  inputted  by  the  user. 


(normally 
set  to  1) 


Standard  printout  (this  includes  a 
call  to  OUTPUT  after  the  solution  of 
every  subproblem) .  Also  the  esti¬ 
mates  of  the  "Lagrange  multipliers" 
and  first  and  second  order  solution 
estimates  are  printed. 

For  additional  printout  (includes 
standard  printout  and  every  interme¬ 
diate  point,  gradient  of  P  and  the 
vector  S) . 


(normally 
set  to  1) 


Final  convergence  is  determined  on 
the  basis  of  current  solution  to  the 
subproblem. 

Final  convergence  is  determined  on 
the  basis  of  the  first  order  esti*- 
mates.  The  first  order  estimate  of 
the  solution  vector  must  be  close  to 
feasible.  See  below. 

Final  convergence  is  determined  on 
the  basis  of  the  second  order  esti¬ 
mates.  The  second  order  estimate  of 
the  solution  vector  must  be  close  to 
feasible  before  the  convergence 
check  is  made.  If  x  is  a  solution 
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Option 

Column 

Value 

Meaning 

estimate  it  is  considered  close  to 
being  feasible  if  g^(x)  +  9o  >  0, 
i=l,2,...,m,  '..’here  60  is  defined  on 
the  parameter  card. 

5 

35 

The  convergence  criterion  determin- 

(normally 

ing  the  NLP  problem  has  been  solved 

set  to  1) 

(only  use  =1,  when  NT4/1) . 

(see  (10, 
p.  195]) 

=1 

Quit  when  — ,G  .  — rr  <  So 

G[x(rk) ,M(rk) ,X(rk)]  0 

=2 

Quit  when  |r  [?=1  £ngj[x(rk)]|  <  60 

=3 

Quit  when 

first  order  estimate  of  vn 

Gfx(rk),y(rk),A(rk>]  0 

6 

42 

=1 

After  final  convergence  the  program 

(normally 

reads  in  new  data  and  solves  the 

set  to  1) 

=2 

next  problem. 

After  final  convergence  has  been  de¬ 
termined  a  call  to  PUNCH  is  made  be¬ 
fore  proceeding  on  to  the  next  prob¬ 
lem. 

7 

49 

First  move  after  a  minimum  to  a  sub- 

(normally 

problem  is  achieved 

set  to  1) 

=1 

No  extrapolation. 

=2 

Extrapolate  through  last  two  minima. 

=3 

Extrapolate  through  last  three  mini- 

ma. 

8 

56 

Not  used. 

9 

63 

=  1 

Subproblem  convergence  criterion,  or 
when  to  stop  minimizing  P  function 
for  fixed  value  of  r  (see  parameter 
card) . 

Quit  when 

|vX(«V)j\-^p]  Vxw(x>,r)|<  r. 

Quit  when 

lVxWT(x1,r)[^]  Vxw(x\r)| 

=2 

„  W (xi_1 )  -  W(xi) 

5 
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Table  3 — continued 


Option 

Column 

Value 

Meaning 

=3 

Quit  when  |VxW(x*,r) |  <  e 

ioa 

70 

=1 

At  least  one  nonlinear  constraint 

=2 

Linear  constraints 

=  3 

Linear  constraints  and  linear  objec¬ 
tive  function  (i.e.,  a  linear  pro¬ 
gramming  problem) 

^hen  option  10=3,  MATRIX  (the  user  subroutine  supplying  the  second 
partial  derivatives)  will  not  be  called,  and  when  option  10=2  it  will 
be  called  only  to  get  the  second  partials  of  f(x). 


TABLE  4 


TOLERANCE  CARD 


Columns 

Format 

Name 

Description 

01-12 

E12.6 

XEP1 

If  some  first  or  second  derivatives  axe 
to  be  gotten  by  numerical  differencing 
(see  addition  option  card)  this  is  the 
value  used  to  compute  them.  See  the 
description  of  DIFF1.  (Usually  setting 
XEP1  equal  to  .0001  is  satisfactory, 
although  a  good  value  is  dependent  on 
the  scaling  of  the  problem.) 

13-24 

E12.6 

XEP2 

When  minimizing  tiie  w  function  for  a 
given  value  of  r  (RH0)  the  value  of  W 
must  decrease  by  an  amount  exceeding 

XEP2  for  each  iteration  after  the 
first.  If  it  does  not,  then  the  code 
prints  out  the  message  "apparently 
roundoff  errors  prevent  a  more  accurate 
determination  of  the  minimum  of  this 
subproblem,"  and  it  is  assumed  a  mini¬ 
mum  has  been  found.  (Usually  we  set 

XEP2  equal  to  0.) 
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TABLE  5 


SECOND  OPTION  CARD 


Value 


Meaning 


=1 

=2 

=3 

=4 

=5 


Solve  problems  without  checking  de¬ 
rivatives 

Solve  problem  after  checking  only 
first  derivatives 

Do  not  solve  nroblem  after  checking 

only  first  derivatives 

Solve  problem  after  checking  first 

and  second  derivatives 

Do  not  solve  problem  after  checking 

first  and  second  derivatives 


=1 


=2 


=3 


=4 


The  method  for  minimizing  the  uncon¬ 
strained  penalty  function  is  to  be 
the  generalized  Newton-Raphson  meth¬ 
od  as  modified  to  handle  indefinite 
Hessian  matrices.  This  method  re¬ 
quires  function  values,  first  and 
second  derivatives. 

Same  as  1,  except  that  when  an  "or¬ 
thogonal  move"  is  made  because  of  an 
indefinite  Hessian  matrix,  -VP  is 
added  to  the  orthogonal  move  vector. 
Steepest  descent  is  used  to  minimize 
P-function. 

The  method  for  minimizing  the  uncon¬ 
strained  penalty  function  is  McCor¬ 
mick’s  modification  of  the  Fletcher- 
Powell  method  us  reported  in  [10]. 
This  requires  function  values  and 
first  derivatives. 


=0 


Do  not  conduct  a  sensitivity  analy¬ 
sis. 

Conduct  a  sensitivity  analysis  at 
the  final  subproblem  with  a  fixed 
value  for  the  differencing  interval. 
Conduct  a  sensitivity  analysis  at 
each  subproblem  along  the  minimizing 
trajectory  with  the  value  of  the 
differencing  interval  depending  on 
the  particular  subproblem. 

Conduct  a  sensitivity  analysis  at 
the  final  subproblem  for  a  range  of 
differencing  intervals. 


Table  5 --continued 


Value 

=0 

=  1 

=0 

=1 

=2 

=0 

=1 

=1 

=2 

=3 

=4 

=5 

*6 


Meaning 


Do  not  estimate  the  partial  deriva¬ 
tives  of  the  estimates  of  the  La¬ 
grange  multipliers. 

Estimate  the  partial  derivatives  of 
the  estimates  of  the  Lagrange  multi¬ 
pliers  whenever  a  sensitivity  analy¬ 
sis  of  the  solution  point  is  con¬ 
ducted. 

Estimate  the  oartial  derivatives  of 
the  optimal  value  function  and  elim¬ 
inate  those  parameters  which  do  not 
affect  the  optimal  value  functions 
from  subsequent  sensitivity  calcula¬ 
tions  . 

Estimate  the  partial  derivatives  of 
the  optimal  value  function  with  re¬ 
spect  to  all  parameters,  but  contin¬ 
ue  all  subsequent  sensitivity  calcu¬ 
lations  with  respect  to  all  param¬ 
eters. 

Do  not  estimate  the  partial  deriva¬ 
tives  of  the  optimal  value  function 
first.  Conduct  the  sensitivity 
analysis  with  respect  to  all  param¬ 
eters. 

Do  not  transform  the  results. 

The  problem  being  solved,  P(x),  is  a 
convex  equivalent  of  a  geometric 
programming  problem  G(t)  obtained  by 
transformation  =  e~xi.  Thus  back 
transform  the  results  to  t  space. 

f*(e)  is  convex.  Derive  parametric 

upper  and  lower  bounds  on  f*(e). 

f*(e)  is  convex.  Derive  parametric 

upper  bound  on  f*(e). 

f*(e)  is  convex.  Derive  parametric 

lower  bound  on  f*(e). 

f*(e)  is  concave.  Derive  parametric 

upper  and  lower  bounds  on  f*(e). 

f*(e)  is  concave.  Derive  parametric 

upper  bound  on  f*(e). 

f*(e)  is  concave.  Derive  parametric 

lower  bound  on  f*(e). 


Figure  1. — Data  deck  structure  for  SENSUMT 
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In  this  section  we  provide  the  input  liscing  and  output  descrip¬ 
tion  of  two  illustrative  examples  which  are  taken  from  [11]. 

5.1  Example  1 

Consider  the  following  convex  RHS  programming  problem: 

2  2 

minimize  f(x)  =  (x^-4)  +  (x2~2) 

2 

subject  to  g^x)  =  -x.^  +  x2  > 

g2(x)  =  -  x2  >.  z2  . 

It  is  desired  to  solve  and  analyze  this  problem  for  sensitivity 
when  =  0  and  e2  =  -3  '  Moreover»  it  is  desired  to  derive  paramet¬ 
ric  upper  and  lower  bounds  on  the  optimal  value  function  of  this  problem 
when 

(i)  -1  <  <0  while  e2  =  -3  ,  and 

(ii)  -3  <  e2  <  -1  while  =  0  . 

5.1.1  Computer  listing  of  the  code  deck.  1 

Figure  2  shows  the  computer  listing  of  the  code  of  Example  1. 

The  letters  in  circles  categorizing  the  input  data  correspond  to  those 
indicated  in  the  code  deck  structure  depicted  in  Figure  1.  The  format 
of  the  listed  data  is  given  in  Section  4. 

5.1.2  Selected  pages  from  the  computer  output. 

Figure  3  lists  an  annotated  computer  output  for  Example  1.  An 
explanation  of  the  meaning  of  the  corresponding  steps  follows. 
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scwisiasru 

IfZMaaaj 


0RGI 


0 

IMPLICIT  REAL*8(A-H,0-Z> 

COMMON/SEN /PAR (45 ) , OPAR( 45 ) , NPAR, IS  ENS 
COMMON/ABO/LY , LZ , PERI 45 ) 

NPAR»2 
PARI  I >-0, 

PAR(2>— 3.0 
PERU)—). 

PER(2>-2. 

RETURN 

iS^RtiDT  nt"  GE  STOTT 

IMPLICIT  PEAL*8(A-H,(>-Z) 

COMMON/SriARE/X  ( 45  > , DEL ( 45) ,  A  ( 45 , 45 ) .  N ,  M  ,  MN ,  NP  I ,  K'W  I 
COMMON/SEN/P AR<  45) ,DPAR(  45 ) , NPAR , I  SENS 
COMMON/ ABC/LY.L2. PERI  45 ) 

I-IN+I 

'  GO  TO  (20,1,2),! 

20  VAL»(  X(  I.)— 4.  )**2*(  X(2)-2.  )**2 
RETURN 

1  VAL— X(  I  K-*2«*<2  >-PAR(  I ) 

OPTIIDM 

2  VAL— X(  I  )— X(2)-PAR(2) 

RETURN 

HSrouTIWE'  (Tra^TT  ITT) 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON/SHARE/X ( 45 ) , DEL ( 45 ) ,  A ( 45 , 45 ) , N , M  ,VN , NP I , NM I 
COMMON/SEN/ PAR<45 ) , DPAR( 45 ) , NPAR, ISENS 
J-IN+I 
DO  5  I«I ,N 
DELUJ-O.O 
5  CONTINUE 

GO  TO  (20, I ,2) , J 
20  DEL(I >-2.*(X(l)-4.) 

OEL(2)-2.*(X<2)-2.) 

RETURN 

1  DEL( I )«-2.*X( I ) 

DEL(2)-1.C 

RETURN 

2  DEL( I )•— I .0 
DEL(2)—J.O 
RETURN 

submurrThf  matrTxTTnTkT 

IMPLICIT  REAL*8 ( A— H, O-Z ) 

CO  MMON/SHARE/X  ( 45 ) , DEL ( 45  > ,  A  (  45 , 45 ) ,  N,  M  VMN,  NP  I ,  NM  I 

COMMON/SEN/PAR( 45), DPAR( 45), NPAR, ISENS 

L-IN-H 

GO  TO  (20, I ,2) ,L 
20  A( 1,1 >-2.0 
A(l,2>-0. 

A(2,2 )-2,  0 
RETURN 

1  A(  1,1)— 2.0 
A(l,2>-0. 

A(2,2)«0. 

*  RETURN 

2  K-t 
RETURN 


EXEC  FORG6,DSNid«79966l  .BOUNDS45', PROG-MAIN 
//GO.SYSLIB  OD 
//CD 

//  DO  DSN-ONU.FORTLIB.DISP-SHR 
//Go.SYSLIN  CD 

//  DO  OS N-&TEMPUNCH,DISP-<OLD. DELETE) 


I  WSJ  31  •liaMl 

I aiffli 


Job  cord 


User  subroutines 


iiMM  — 

n^rn-nrrrnrrn: 


[o»  ]->nrnr— I 

I'VPRviojrwi 


Figure  2. — Computer  listing  of  the  code  for  Example  1 
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■Annotated  computer  output  for  Example 


»M*  if  VALUt  Jf  X  IS 

i.iiJiM.M  ai  o«uvr20*o  at 

f H  CJNSTrf&IMT  VALUES 
»•  1450M9J-0J  0.  AS 76* 720-0* 


LAGftAUuC  MULTIPLIERS 

0.  / 1666910  01  P»  0.7)667010  01  C»  0.73666930  01  RSIGXA-  0.0 

!►*  CUKXE'H  value  OF  X  IS 


Figure  3. — continued 
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Figure  3. — continued 
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Figure  3. — continued 
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M 

RESIGMA  =  -r  Jt  £n  g .  [ x  (r>  ]  (note  that  this  may  be  negative) 

j=l  3 

M+MZ  „ 

H  =  l  h^[x(r)]/r 
j=M+l  3 

G  =  dual  value  =  f[x(r)]  -  r*M  +  2.*H 

(in  a  convex  problem  F  _>  v*  >.  G) 

where  the  current  value  of  x  is  x(r)  and  the  current  con¬ 
straint  values  are  g1[(r)]  and  g£[x(r)]  . 


Values  of  F,  W  (designated  by  P  in  the  computer  output),  G, 
RSIGMA,  and  x  are  repeated.  Current  constraint  values  (i.e., 
Uj  and  Wj)  are  r/g^  »  >  and  2n4(x)/r  ,  j=M+l,M+MZ. 

Note:  Since  Example  1  has  no  equality  constraint,  the  entries 
relating  to  h^  in  8  and  9  are  zero. 

The  solution  to  the  second  subproblem  (r  -  10/4  =  2.5). 

The  current  first  order  estimates  of  x  (obtained  by  first 
order  estrapolation,  using  the  solution  of  the  first  and  second 
subproblems)  and  the  corresponding  problem  functions. 


The  solution  to  the  third  subproblem  (r  =  2.5/4  =  .625). 


The  current  second  order  estimates  of  x  (obtained  by  second 
order  extrapolation  using  solutions  of  the  first,  second,  and 
third  subproblems)  and  the  corresponding  problem  functions. 


*1 

i 


The  solution  to  the 
The  solution  to  the 
The  solution  to  the 
The  solution  to  the 
The  solution  to  the 
The  solution  to  the 
The  solution  to  the 


fourth  subproblem, 
fifth  subproblem, 
sixth  subproblem, 
seventh  subproblem, 
eighth  subproblem, 
ninth  subproblem, 
tenth  subproblem. 
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The  solution  to  the  eleventh  subproblem. 

The  solution  to  the  twelvth  (final)  subproblem. 

The  general  information  about  the  sensitivity  analysis  including 
the  value  of  r  ,  the  estimated  final  solution  point  and  the 
parameter  values  and  associated  differencing  intervals. 

Estimates  of  the  sensitivity  of  the  optimal  value  function  ob¬ 
tained  by  taking  the  gradient  of  the  Lagrangian  with  respect  to 
the  parameters  as  described  in  Step  9  of  Algorithm  2.1. 

The  parameters  for  which  detailed  sensitivity  results  follow. 
Here  it  is  both  parameters  since  the  optimal  value  function  is 
sensitive  to  both. 

Estimates  of  the  partial  derivatives  of  the  solution  point  taken 
with  respect  to  parameter  one  as  in  Steps  1-3  of  Algorithm  2.1. 

Estimates  of  the  partial  derivatives  of  the  Lagrange  multipliers 
taken  with  respect  to  parameter  one  as  in  Step  5  of  Algorithm 


Estimate  of  the  partial  derivative  of  the  optimal  value  function 
taken  with  respect  to  parameter  one  and  obtained  by  using  the 
chain  rule  as  in  Step  8  of  Algorithm  2.1. 


Same  as  (26)  ,  (27)  ,  and  (28/  but  with  respect  to  parameter  two. 


Output  for  solution  of  the  first  perturbed  problem,  i.e.,  r.^  = 

-1  ,  e2  *  (Step  2  of  Algorithm  3.1). 

Printout  of  the  input  data  and  subproblem  solutions  for  the 
first  perturbed  problem. 

Printout  of  the  final  subproblem  for  the  first  perturbed  problem. 
Same  as  (2$)  but  for  the  first  perturbed  problem. 


Same  as  (24)  but  for  the  first  perburbed  problem. 
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Output  for  optimal  value  function  bounds  as  a  function  of  first 
parameter  (Steps  3  and  4  of  Algorithm  3.?). 

(f*(e^),e^)  of  the  unperturbed  and  first  perturbed  problems 


[see  (2jj)  ,  (23)  ,  and  (3p  ,  Q2)  ,  respectively]. 

Equations  of  parametric  upper  bound  on  f*(e^)  . 

Equations  of  parametric  lower  bounds  or.  f*(e^)  . 

Maximum  of  these  equations  over  the  range  of  -1  £  <  0  pro¬ 

vides  a  lower  bound  on  f*(e^)  . 

Equation  of  a  quadratic  estimate  of  f*(e^)  [see  the  descrip¬ 
tion  of  subroutine  BOUND  in  Section  6  fcr  the  method  of  deriving 
this  equation] . 

Value  of  bounds  and  the  quadratic  estimate  of  f*(e^)  at  eleven 
equidistant  points  over  the  perturbation  range  of  the  first 
parameter  . 


Output  for  solution  of  the  second  pertuibed  problem,  i.e.,  = 

0  ,  e2  =  -1  (Step  2  of  Algorithm  3.1). 

Description  of  the  output  for  the  second  perturbed  problem  is 
similar  to  that  of  the  first  perturbed  problem,  i.e.,  ©  -  © 


Output  for  optimal  value  function  bounds  as  a  function  of  the 
second  parameter  (Steps  3  and  4  of  Algorithm  3.1). 

The  description  here  again  closely  parallels  that  of  the  first 
parameter,  i.e.,  ©  -  ©• 


The  graphical  depiction  of  the  bounds  derived  for  f*(e)  as  a 
function  of  »  together  with  a  plot  of  the  analytical  solution  and 

quadratic  estimates  of  f*(e)  at  different  value.t  of  is  in  Figure 
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Figure  4. — Graph  of  bounds  on  f*(e2>  of  Example  1  (computer  solution) 
5.2  Example  2 


Minimize  f(x)  =  x2  +  .5e  2  -  ,125x^  +  ,25(x  -40)2  -  .le  5 

» • • . , Xj 

Subject  to  g1(x)  -  -.5x2  +  6x2  -  5e  3  -  .05x2  -  .5/x^  > 

~X1  2 

g2(x)  *  -5e  -  2x2  +  3x3  +  y^  +  3x5  >  -12 


g3(x) 


2  2 

3x.  +  x„  -  x,  +  r. ,  -  xc  >  2 

12  3  ‘t  5  — 


0  <  x.  <  5  ,  0  <  x,  <  5  ,  0  <  x  <  •  ,  j-1,2,4 
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It  is  desired  to  calculate  bounds  on  f*(e^)  of  this  problem  for 
eL  e  [-10,5]  . 

(i)  Problem  SR(z). 

•  Standard  format: 

5 

Min  f(x)  =  [  f,(x.) 

xlt...,x5  j*l  3  3 

S.t.  x  e  G  H  B  ,  where 

G  =  (x  :  ,  g2(x)>-12  ,  g^(x)>2}  ,  -10<e1<5 

and 

B  =  (0£x^<^5  ,  0£x^<^5  ,  0  <  x_.  5  °°  ,  i*l,2,4)  . 

(ii)  Problem  CSH(e),  (convex  overestimating  problem  of  SR(z.)). 

•  Problem  formulation: 
f3(x3)  =  -3.125x3  +  6.014 

f5(x5)  =  -2.948x5  +  7.027 

Remaining  terms  are  identical  to  those  of  problem  SR(e); 
thus  problem  CSR(e)  reads 


2  2  3 

*  x^  +  .5e  -  .125x3 


2  5 

+  .25(xa-40)  -  .le 


Min 

» *  *  *  i  ^ 


f(x)  =  x2  +  .5e  2  -  3.125x3  +  .25(xA-40)2 
-  2.948x5  13.041 


s.t.  x  e  5  n  b  ,  -lo^e^s  , 

where  8  =  G  . 

Note:  f3(x3>  and  f^(x^)  were  chosen  here  to  be  the 

lowest  nonunderestimating  lines  parallel  to  the  convex 
envelopes  of  f ^ (x^>  and  f3(Xj)  »  respectively. 
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(Hi)  Problem  CSR (c) ,  (convex  underestimating  problem). 

•  Problem  formulation: 

=  -3.125xj 

f^(x,.)  =  -2.948x,. 

Remaining  terms  are  identical  to  those  of  problem  SR(c); 
thus  problem  CSR(e)  reads 

Min  f(x)  =  x2  +  .5e  2  -  3.12x3  +  .25(x4-40)2 

» • • • * 

-  2.948x5  -  .1 

S.t.  x  e  G  O  B  ,  -10<  £x<  5 

where  2  =  G  . 

The  computer  listing  of  the  code  for  problem  CSR(e)  is 
shown  in  Figure  5.  The  listing  of  the  code  for  CSR(e) 
is  identical  to  that  of  the  problem  CSR(e)  except  for 
the  constant  term  in  the  objective  function,  which  is 
+  13.041  ,  rather  than  -  .1  . 

(iv)  Bounds. 

Figures  6  and  7  depict  the  calculated  upper  and  lower 
bounds  via  the  analysis  of  the  problems  CSR(e)  and 
CSR(e),  respectively.  These  results  are  summarized 
below. 

Upper  bound: 

f*^)  =  4.825e1  +  136.628  . 

Lower  bound: 

f*^)  *  max[1.880ei  +  94.032  ,  +  120.284]  . 

The  graphical  depiction  of  these  bounds  is  in  Figure  8. 
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J«GmA£MI  T-2  L*2 

//  EXEC  FORGJ 

SU0WOUT  I  NE  HEAD  IN 
IMPLICIT  REAL*8(A-H.O-Z) 

COMMON/SEN/PAR (  45) ,  DPAR(  45) ,  k(PAR,  ISENS 
GOMMON/ABG/LY, LZ , PER( 45 ) 

NP AR=3 
PARC  I  )*5. 

PAR(2  )=- 1 2 . 

PAR(3)=2. 

PE3( I )=- 1 5. 

RETURN 

END 

SUBROUTINE  RESTNT( IN.VAL) 

IMPLICIT  REAL*8( A-H.O-Z ) 

CD  MMuN/SHARE/X ( 45 ) , DEL ( 45 ) ,  A ( 45 , 45 ) , N, M ,«N, NPI , MM  I 
COMMON/SEN/P AR <  45 ) ,  QnAR  (  45) ,  NPAR,  I  SEN'S 
C0MM0N/ABG/LY.LZ.PER(45) 

1  =  I N+ 1 

GO  TO  (20,1 ,2, 3, 4, 5), I 

20  VAL*X ( I )**2* . 5*DEXP(X(?))-3. I25*X( 3 )*.2b*X X( 4)-40.  >**2 

★  -2.948*X(5)-H  3.041 
RETURN 

1  VAL*-.5*X( I )**2+6.*X(2)-5.*DEXP(X(3>>-.05*X(4)+*2-.5/X(5) 

*  -PAfl( I ) 

RETURN 

2  VAL— 5.*DEXP(-X( I ) )-2.*X ( 2 ) **2+3 .*X ( 3 ' +X (4 )+3.*X ( 5 )-»AR( 2 ) 
RETURN 

3  VAL-3.*X( I )*X(2 )-X( 3>**2*X(4)-X(5>**2-PAR(3> 

RE  riJPN 

4  VAL*- X(3)*5. 

RETURN 

5  VAL»-X(5)+5. 

RETURN 

END 

SUBROUTINE  GPADI ( IN) 

IMPLICIT  REAL*8(A-H,0-Z> 

COMMON/SHARE/X  (  45 ) ,  DEL  (  45 ) ,  A  ( 45 , 45 ) ,  N,  N  ,MK ,  NPI  ,  NM  t 
CoMMON/SEN/PAR( 45) , DPAP<  45) , MPAP, ISENS 
J=  I  N«-  I 
DO  15  I*I,N 
DEL< I )-0.0 
15  CONTINUE 

GO  TO  (20, 1,2, 3,4,5), J 
OEL (l>-2.*X(I) 

('EL (2  )«.  5*DEXP(  X(  2  ) ) 

DEL(3)»-3. (25 
DEL(4)=.5*(X(4)-40.) 

DEL<5)«-2.948 

RETURN 

1  l)EL(l  )  — l.*X(  I  ) 
l)FL(2>*6. 

i)EL(3)*-5.*I)EXP(  X(  3 )) 

DEL<  4 )*— . I *X ( 4 ) 

OEL( 5 )  =  . 5/X<  5) **? 

RE  TURN 

2  )EL( I )=5.*DEXP(-X( I >) 

DE'.(2)=-4.*X(2) 

>EL(  3 )  =  3 . 
i'-:r.(4)*i, 

5  )  =  3. 


Figure  5. — Computer  listing  of  the  code  for  the  convex 
overestimating  problem  of  Example  2. 
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re  rur.fi 

1  0ELU>  =  3. 

OEL(  2  )  =  l  . 

DEL(3>— 2.*X(3> 

DELC4 )*1 . 

DEL(5)*-2.*X(5> 

RETURN 

4  DEL<  3 )  =  - 1 . 

RE  TURN 

5  0EI.(5)=-I . 

RETURN 

END 

SUdROUTINE  MATRIX(IN.K) 

IMPLICIT  REAL*8 ( A-H, O-Z > 

CO'(WON/SHARE/X  <  45 ) «  DEL ( 45 )  ,  A < 45 , 45 > . N, M ,MN, NPI  ,N'<| 
COMMON/SEN/PAR<  45  > .  DPAR  ( 45 ) ,  NPAR,  ISE'JS 
L* IN+ I 

GO  TO  (20, 1,2, 3, 4, 4), L 
20  A<  1,1 ) =2 . 

A(2,2)=.5*DEXP(X(2)> 

A(4,4)=.5 

RETURN 

1  V(l, )>*-!. 

A( 3, 3 )=-5.*DEXP( X( 3 ) ) 

A<  4,4 )=-. I 

A( 5,5 )*-l . /X ( 5 ) **3 
RETURN 

2  A( 1,1 )=-5.*DEXP(-X( I)) 

A(2.2)*-4. 

HE  TURN 

3  A(3,3)=»- 2. 

A( 5,5 I=-2. 

RETURN 

4  K*l 
RETURN 
LNO 

//  EXEC  EoRG6,0SN--'OR7995<51  .B0UNDS45' . PROG*MAI\' 

//Go.bYbLIB  00 

//  no 

//  00  05N=GrtU.FoHTLI6,DISr’=SHR 

//GO.SYGLIN  DO 

//  00  !)5N=4TEMPUNCH, 0 I SP=( OLD, DELETE) 

//F107FO0  I  no  SYSfXJT=A 

I.0E-08  10. 0  I.0E-06  4.0  900.0  5 

3.  3.  4.  I.  4. 

3  I  I  I  I  I  )  1 

0. IE-08  0.0001 

14)1102 

// 
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UPIMAL  VALUE  F  UNt,  I  I  UN  aUJNUS  WHEN  PAKI  1)  IS  PERTURBED 

*v*******t*«*»tt ***** **»***<.***»«»*»***** **»»****»***•»•• 


Plllll  l  (UNPLR  TURDE  J  SOLUT  I  UN  I  : 

PAW  I  Ih  O.SOOUOOO)  01  T  (PAM  III-  0.160  754  71)  03 


PDIi. I  2  IPEKTU.OiEJ  SULUTIUN)  : 

p  U.  (  ni--  -).  I  MOODOI)  02  F  t  P  \R  I  11=  a.fl837)77U  02 


till  Prt'.ilOO  TilPUlir.ll  PUINIS  1  A1)')  2  AN')  UVIR  ESTI  M.’.  T  IN  i  f  * 


>.  *125  3000  01  *  >Akl  11  ♦  0.13662770  0  3 


l<  ,,,))!)  tVALOUI'il  At  TEN  f.)UI  al  S  f  ANCF  POINTS  BETWEEN  POINTS  1  ANO  2 


►Mi  (  l)  LOWER  UUUNO 


UPPER  BOUND 


QUAD  ESTIMATE 


4.0)0 
».  >00 
••'I  )0 
I.  400 
-1.000 
-2.4  1  ) 
-  *.000 
- >.  500 
-  7.000 
-I.  6  JO 
- 1*.|)0  * 


0. 16075470  03 
0.  1 5 35 ’.661)  03 
0.1462/R50  03 
0. 1 300**341)  03 
0.13160230  03 
0. 12457430  03 
0.11732621)  03 
0. 1 1006810  03 
0.1 023 5001)  03 
0.9561  >6  70  02 
0.68373770  02 


Figure  6. — Parametric  upper  bound  on  f*(E^) 
Example  2  (computer  solution). 


via  problem  CSR(e), 
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(Vli;ui  VAl'JI  FUNCTION  I30UNOS  *HLN  PARI  II  IS  PERT'JRflEO 


euriT  I  I  JUPLul ’J:UU  0  SULUTlUitl  : 

PAUI  11=  v1.  5  UUOOUJJ  0 1  I  I  PAR  f  II  )=  0.t<t7ol3T0  03 

Pill'll  2  (PI.KTUKMF)  SOLUTION!  : 

PAi  I  111=  -O.iOOJOOOU  02  f  I  PAR  I  11=  0.75232:70  02 


tin.  LSIlilAIMo  I*  AT  PUl  If  l 


I  =  0.  V.6SV2I  i)  01  *  PAR  I  II  ♦  0.  12  02  A  All)  03 


MOt  0101k  iSIMAHUii  1*  A 1  PU I  NT  2 


0.1«l')9StU  H  ♦  PARI  n  ♦  0.9A032SRI3  02 


ft  ,  HID  I  VALUATION  AT  ILN  F'JU  1 01S  T  4NCE  POINTS  BETWEEN  POINTS  1  ANO  2 

PAR (  II  LOWER.  HOUND  UPPER  ROUND  OUAO  ESTIMATE 


3.  >00 

o.i-»ri>ii?o 

03 

300 

O.H9'<IAOO 

01 

2.0  00 

0.  1  11 21. SOU 

03 

J.  ‘.o' I 

0.  12  301  no 

03 

-I. 00) 

0.  1146 1020 

01 

-2.30 ) 

0. 106010  JO 

03 

J  JO 

0.91A2JA60 

02 

-■».  300 

0.0022  !'•  80 

02 

-  1. 000 

0.R2J2270U 

02 

-H. 500 

J.  7H03275D 

02 

-10.000 

0.  15212  no 

02 

Figure 
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Parametric  lower  bound  on  f*(e^) 
Example  2  (computer  solution). 


via  problem  CSR(c), 


\ 
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Figure  8.  —  Parametric  bounds  on  f*(e)  ,  Example  2 
(computer  solution). 
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6.  General  Description  and  Listing  of  the 
SENSUMT  Subroutines 


The  subroutines  comprising  SENSUMT  fall  into  four  categories: 

(i)  user  subroutines, 

(ii)  SUMT  subroutines, 

(iii)  sensitivity  subroutines,  and 

(iv)  bound  subroutines. 

As  mentioned  before,  user  subroutines  are  generally  composed  of 
four  subroutines,  i.e.,  READIN,  RSTNT,  GRAD1,  and  MATRIX.  These  subrou¬ 
tines  respectively  provide  pertinent  information  about  the  problem, 
functions  of  the  problem,  their  gradients  and  Hessian  matrices.  A  de¬ 
tailed  description  of  these  subroutines,  together  with  the  corresponding 
codes  for  Examples  1  and  2,  were  given  in  Section  4  and  Figures  2  and  5. 

SUMT  subroutines  implement  the  Sequential  Unconstrained  Minimiza¬ 
tion  Technique  of  Fiacco  and  McCormick  [10].  The  subroutines,  along  with 
program  MAIN,  comprising  this  group  are: 


•  B0DY 

• 

0UTPUT 

•  CHCKER 

• 

PEVALU 

•  C0NVRG 

• 

PUNCH 

•  DIFF1 

• 

REJECT 

•  DIFF2 

• 

RH0C0M 

•  ESTIM 

• 

SEC0ND 

•  EVALU 

• 

SEC0RD 

•  FEAS 

• 

SET 

•  FINAL 

• 

ST0RE 

•  GRAD 

• 

TCHECK 

•  INVERS 

• 

TIMEC 

•  MAIN 

• 

TRANS 

•  0PT 

• 

XM0VE 

With  the  exception  of  the  subroutine  TRANS,  the  above  routines 
have  been  developed  by  Mylander,  Holmes,  and  McCormick  [13].  Some  of 


-  59  - 


T-434 


* 


these  routines,  in  particular  program  MAIN  and  subroutines  B0DY  and 
0UTPUT,  have  been  modified  for  implementing  sensitivity  and  bound  calcu¬ 
lation  routines.  Subroutine  TRANS  was  recently  developed  to  aid  the 
user  when  analyzing  a  convex  equivalent  of  geometric  programming  prob¬ 
lems  by  SENSUMT.  See  page  104  for  a  more  detailed  description  of  this 
subroutine. 

Sensitivity  subroutines  implement  the  sensitivity  analysis  Algo¬ 
rithm  2.1  and  interface  it  with  the  SUMT  subroutines.  A  brief  history 
of  the  development  of  these  subroutines  was  given  in  Section  1.  The 
latest  version  of  the  codes  comprising  this  group,  due  to  Armacost  [1], 
are  subroutines  SENS,  PARDIF,  LMULT,  and  PRESEN.  Subroutines  SENS  and 
PRESEN  have  been  slightly  modified  in  implementing  the  bound  calculation 
routines . 

Bound  subroutines,  developed  in  [11],  implement  the  bound  calcu¬ 
lation  Algorithms  3.1  and  3.2  given  in  Section  3.  The  two  subroutines 
in  this  group  are  BOUND  and  PERT. 

All  of  the  subroutines  comprising  SENSUMT  are  dimensioned  to 
solve  problems  having  at  most  45  variables,  45  parameters,  and  200  con¬ 
straints.  However,  if  the  computer  capacity  permits,  they  may  readily 
be  re-dimensioned  to  solve  problems  of  larger  size.  All  of  these  sub¬ 
routines  are  separately  filed  in  the  Conversational  Monitor  System  (CMS) 
component  of  the  IBM  Virtual  Machine  facility  37C  (VM/370)  at  The  George 
Washington  University  Center  for  Academic  and  Administrative  Computing 
under  their  corresponding  names. 

To  make  this  manual  self-contained,  the  computer  listing  and  a 
general  description  of  each  subroutine  is  provided  (in  alphabetical 
order  by  name).  The  descriptions  of  the  SUMT  subroutines  are  largely 
taken  from  [13],  and  those  of  the  sensitivity  subroutines  are  taken 
from  [1].  The  user  subroutines  are  problem-dependent,  thus  they  are  in¬ 
cluded  in  the  following  listings.  For  familiarity  with  the  user  subrou¬ 
tines,  the  reader  should  refer  to  the  coding  of  Examples  1  and  2  pro¬ 
vided  in  Figures  2  and  5,  respectively. 
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6.1  B0DY 

Subroutine  B0DY  coordinates  the  flow  among  the  subroutines  that 
actually  do  the  calculations  required  by  the  various  phases  of  the  al¬ 
gorithm.  The  flow  in  this  routine  is  slightly  different  when  the  pro¬ 
gram  is  in  the  feasibility  phase  (solving  the  entry  problem)  rather  than 
the  normal  phase  (solving  the  stated  NLP  problem) .  The  listing  appears 
as  Figure  9. 

6.2  B0UND 

The  subroutine  B0UND,  called  by  the  program  MAIN,  generates  the 
equations  for  the  upper  and  lower  bounds  of  the  optimal  value  function 
as  functions  of  the  problem  parameters.  When  analyzing  a  problem  which 
is  known  to  have  a  convex  or  concave  optimal  value  function,  B0UND  also 
derives  the  equation  for  a  quadratic  function  which  approximates  the 
optimal  value  function.  If  f*(e)  is  convex  (concave),  it  corresponds 
to  the  higher  (lower)  of  the  two  quadratic  functions,  in  the  perturba¬ 
tion  range  of  the  parameter  of  interest,  which  passes  through  the  two 
calculated  points  in  f*(e),e  space  and  has  the  calculated  slope  at 
each  of  these  points.  Evaluation  and  printout  of  this  quadratic  func¬ 
tion  (when  applicable),  and  the  upper  and/or  lower  bound  at  eleven  equi¬ 
distant  points  over  the  perturbation  range  of  the  parameters  of  interest 
are  also  programmed  in  subroutine  B0UND.  The  listing  appears  as  Figure 
10. 


6.3  CHCKER 

Subroutine  CHCKER  evaluates  the  first  partial  derivatives  for  all 
the  functions  at  the  starting  point  first  by  calling  the  user-supplied 
subroutine  GRAD1  and  then  by  calling  DIFF1.  The  results  are  printed  out 
to  aid  the  user  in  debugging  the  subroutines  used  to  describe  the  NLP 
problem  he  wishes  to  solve. 

The  matrix  of  second  partial  derivatives  of  each  function  is  also 
evaluated  by  the  two  methods;  first  by  calling  MATRIX  and  then  calling 
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0031 

SUBROUTINE  BODY 

§0090363 

C 

8003037 ) 

c 

AUGUST  1971 

§09000*0 

c 

§0000090 

:  BODY  COORDINATES  THE  FISH  AMONG  THE  SUBROUTINES  THAT  ACTUAL  IV  DO 

THE  80903100 

C  CALCULATIONS  RE  OU 1 R  E  0  BY  THF  VARIOUS  FARTS 

OF  THE  ALGORITHM. 

80930119 

000? 

IMPLICIT  AEALFBIA-H.0-21 

80090129 

003  4 

REAL *4  HHOt  N , RATIO, EFSI.T HE TAO.Af FI »  XF  P2 

800091 33 

00  )<. 

C0MM0N/SHARE/XI45),0ELI45I,AI45,'.SI  ,N,M 

lfNN,NPl,Nm 

80009140 

000  s 

COMMON  /CPTNS/  NT  1, NT 2, NT3, NT  4, NT  5, NT  6, 

NT  7«  NT8#NT9*NT10 

80000159 

0006 

COMMON/V AL UE/E, G.PO, RSI GMA.AJ 1901 ,RNO 

§0000163 

000 1 

C3MM3N/CRST/0EL  X 1 45 1 ,OELXOf  4  5 1 ,RHOI N, RATIO, EP SI, THE T AO, 

§0000170 

IRSIG1,G1,X1I45I ,X2 1451, X 31 451  , XR21 4  51, XRI 1451, FA  1, 

800001*3 

2FR2,P|,F1,RJ II 901,001 T.FGRAOI 451,01  AG  14  61 , 

§0000190 

3  FKEV3.A0ELX,  NTCTR,  NUNINI,  NPHASE,  NSAT IS 

§0000200 

00  oh 

C0HM3N  /CONRAR/  Nfl*NF2*NF3 

§0000210 

000^ 

CONMON/EX>>OPT/NEXOFL,NE<OF2,NFXCF3.NEXOF4,NEXOFS,XEF1,XEF2 

80000229 

00 10 

CONMON/SEN/PAR  1451,  DPAR  1 45  1 ,  NFAR,  1  SENS 

80090239 

001 1 

COHLQN/ABG/LYtLZ , PER 145 1 

80000249 

OUi? 

C0MMIJN/ABGL/FE1.FE2 

80000259 

0013 

C0MMUN/ABG2/0F 1 1 451 ,DF2I 451 

80090269 

0014 

NE2*2 

80000273 

COl* 

NF3*2 

80000280 

0016 

HN*0 

80000299 

0017 

NUNINI -0 

80000100 

001 * 

1  SENS  •  0 

80000310 

c 

OPTION  OF  GETTING  INITIAL  RHO 

80000320 

0019 

CALL  RHOCOH 

B0000339 

0320 

CALL  EVALU 

§0000140 

0021 

10 

CALL  XNOVE 

§0000350 

002? 

GO  TO  130,201,  NT 3 

§0900360 

0023 

23 

CALL  TINEC 

80000370 

0024 

CALL  OUTPUT  111 

800003*0 

002b 

GO  TO  40 

80000390 

0026 

30 

CALL  TCNECK 

80000400 

C  II 

N  FEASIBILITY  PHASE  4  MEANS  *EAS  ACHIEVED 

80000410 

002  7 

40 

GO  TO  ISO, 50, 50, 2001,  NSATIS 

80000420 

002* 

50 

CALL  CONVRG  INI  1 

80000430 

002* 

GO  TO  160,10.1251,  N1 

80900440 

c 

MINIMUM  ACHIEVED  IF  Nl»l 

80000450 

003  ) 

60 

GO  TO  I/O, SOI,  NTS 

80000460 

0031 

70 

CALL  TIMEC 

80000470 

0032 

CALL  OUTPUT  111 

80000400 

c  - 

—  NUMBER  OF  NINIMA  ACHIEVED  INCREASED  BY 

1 

80000490 

OOii 

40 

NUM1NI>NUN!NI*1 

80000500 

0034 

HN*3 

80000510 

003b 

GO  TO  1190,90,901,  NPMASE 

80000520 

0036 

90 

CALL  ESTIM 

§0000530 

C 

FINAL  MIGHT  HAVE  BEEN  CALLED  BT  ESTIM-CONVERGEO  IF  N2*l 

80000540 

003  7 

IEINEX0P3.NE.2I  GO  TO  27 

80000550 

003b 

CALL  SENS 

80000560 

003v 

27  GO  TO  1100,110,1201.  NT4 

80000570 

C 

NT 4« l  FINAL  CONVERGENCE  ON  0  ORDER  ESTIMATES,  NT4-2  CONVERGE  ON 

FI8SB0D005B0 

C 

OPOER  ESTIMATES,  NT4>3  CONVERGE  ON  SECOND 

ORDER  ESTIMATES. 

80000590 

0060 

100 

CALL  FINAL  INFli 

80000609 

0041 

GO  TO  1 1 30, 1A01 ,  NF1 

80000610 

0042 

110 

GO  TO  1130,1401,  NF2 

80000620 

0043 

120 

GO  TO  1130,1401,  NFI 

80000630 

0044 

125 

NPHASE-5 

80000640 

004  5 

130 

RETURN 

80000650 

0046 

140 

RH0*RH0/RATI0 

80000660 

C  USING  PREVIOSLY  CONFUTED  VALUES  FOR  F,  ANO 

RJ  P  IS  RECONPUTEO  WITH  THER000067D 

C  NEW  VALUE  OF  RHO. 

BOD006R0 

0047 

CALL  PEVALU 

60000643 

C  A 

VECTOR  IS  LEFT  IN  DELXIII  BV  ESTIM 

60000700 

004* 

IT  1NUM1NI-2)  10,150,150 

B0D007I9 

0049 

150 

GO  TO  110,160,1601,  MTT 

BOD00720 

0050 

160 

CALL  GRAO  121 

60000730 

0051 

CALL  OPT 

60000740 

005? 

GO  TO  1  ISO, 1 701 ,  NTS 

60000750 

0053 

l  TO 

WRITE  16,2101 

B 0000 760 

0054 

CALL  OUTPUT  III 

B 0000770 

005S 

ISO 

GO  TO  50 

§0000760 

0056 

1 40 

IF  IGI  90,90,200 

60000790 

005  7 

200 

RETURN 

BOOOOBOO 

C  -  DUAL  VALUE  GREATER  THAN  0  MEANS  NO  FEASIBLE  FOINT  EXISTS 

BOOOOSIO 

c 

60000R20 

0056 

210 

FORMAT  I6X, 30HN0VE0  ON  EXTRAPOLATION  VECTOR  1 

B 0000630 

005'* 

€N9 

B0000640 

Figure  9. — Subroutine  B0DY. 
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C 

SUBROUTINE  BOUND  HAS  DEVELOPED  BY  GHAENII1979I  TO  CALCULATE 

60U00060 

c 

UPPER  ANO/OR  LOVER  P1ECEH1SE  LINEAR  PARAMETRIC  BOUNDS  ON 

BQU00070 

c 

DPI  INAL  VALUE  FUNCTIONS  HHEN  THESE  FUNCTIONS  ARE  KNOHN  TO  BE 

SOUOOO0O 

c 

CONVEX  UR  CONCAVE  OVER  THE  RANGE  OF  A  GIVEN  PARAMETER. 

BOJ 00090 

0001 

SUBROUTINE  BOUND 

60000100 

0002 

IMPLICIT  RE ALPS 1 A-H. 0-2 1 

60U001 10 

0003 

COMMON /SEN/ PARI 45 1 , OP AR (45 1 »NPAR. 1SENS 

60000 120 

0004 

COMMON /ABG/LY .L2.PER145I 

60U00133 

0005 

C0MNUN/ABG1/FE1.FE2 

0OUOO14O 

0006 

C0MMUN/ABG2/0F 1 ( 451 ,DF2 1 451 

SOUOOl SO 

0007 

COMMON/OP7/NEXOP7 

B0U0016J 

0008 

DIMENSION  E (201 .FQI20I ,FUQI 201, FU12 1201 .FL12I201 

Buoooiro 

0009 

El-PARILZi 

60U00160 

0010 

E2-PARIL2 1  *  PER1LZI 

6 OU 00190 

c 

LINE  FIT  TROUGH  POINTS  1  ANO  I 

B0U00200 

0011 

S12  -  1FE1-FE2I/IE1-E2I 

60U00210 

0012 

B12  •  FEI-S12PE1 

BOU00220 

c 

TANGENT  LINE  UNDER  ESTIMATING  FP  AT  POINT  1 

60U00210 

0013 

S1-0FML2I 

0OUOO24O 

0014 

B1»FE1-SI*E1 

00000250 

c 

TANGENT  LINE  UNDER  ESTIMATING  FP  AT  POINT  2 

BOU00260 

0015 

S2-DF21L2I 

0OUOO27O 

0016 

B2-FE2-S2PE2 

B0U00280 

c 

INTERSECTION  OF  THE  ABOVE  TUO  LINES 

BOU00290 

0017 

E3-(B2-S1I/(S1-S2I 

B0U00300 

0018 

FE3-SIPIB2-B11/ISI-S2I  ♦  B1 

80U00310 

c 

QUADRATIC  FIT  TROUGH  POINTS  1  ANO  2 

B0U00320 

0019 

A»(FE2-E2*Sl-f El*El*Sll /( (El-E2l*»2l 

60000330 

0020 

C-((ElP*2IP(FE2-£2PSllP|FEl-ElPSl»P(E2PP2-2.PElPE2l!/( IE1-E2IPP2I 

B 0000 340 

0021 

B-SI-2.PEIPA 

B0U00353 

0022 

AA-( FE1-£1*S2-FE2pE2*S2 1/1 I E2-E1 |P»2I 

B0U0036D 

0023 

CC><IE2P*2|P|FE1-E1PS2)P(FE2-E?PS2)P|E1PP2-2.*E2PE1)1/IIE2-E1IPP2IBOU00370 

0024 

BB«S2-2.*E2*AA 

BOU00380 

c 

BOUND  EVALUATION  AT  10  DISCRETE  POINTS 

B0U00393 

0025 

El 11-El 

B0U00400 

0026 

0E>PER|LZI/10. 

B0U0041 0 

0027 

00  580  1-1,10 

BOU00420 

0028 

IFll.EO.il  GO  TO  750 

B0UD0430 

0029 

Elll-Ell-llOE 

B0UD0440 

c 

QUADRATIC  ESTIMATE 

BOU00450 

0030 

750  FQll l«AP£t !>p*2»BPE1II*C 

BOU 00460 

0031 

FQOI I 1 -AAPE ( I 1 PP24BBPEI 1 1 PCC 

BOU00470 

c 

UPPER  LINEAR  BOUND 

BOU00480 

0032 

FUL21 1  1-S12PEI  1 14BV2 

BOU00490 

c 

LOVER  LINEAR  SOUND 

BOUOOSOO 

0033 

IF1PERIL2I.GT.0.I  GO  TO  5B5 

BOU00510 

0034 

IFIEl II.LT.E3I  GO  TO  570 

80U00520 

0035 

FL12III-SIPEI1I4B1 

BOU00530 

0036 

GO  TO  595 

BDU00540 

0037 

570  FL12I1I«S2PEIII»B2 

B0UD0550 

0038 

GO  TO  595 

B0U0O560 

0039 

585  IFIEl 11. GT.E3I  GO  TO  590 

B0U00570 

0040 

FL1211  I-S1PEII1PB1 

BOUOOSBO 

0041 

GO  TO  595 

80U00590 

0042 

590  FL12III>S2PEI1I*B2 

BOU00600 

0043 

595  CONTINUE 

B0U00610 

0044 

580  CONTINUE 

SOU 00620 

c 

CHOICE  ON  PROPER  QUADRATIC  FIT 

BOU 00430 

0045 

IFINEX0P7  .GE.41  GO  TO  810 

80U00640 

0046 

IFIFQI5)  .GT.  FQQI5II  GO  TO  BOO 

80U00650 

0047 

A-AA 

BOU00660 

0048 

B-BB 

B0U0O670 

0049 

c-cc 

B0U00680 

0050 

DO  100  1-1,10 

80U00690 

0051 

FQIII-FQQIII 

80U00700 

0052 

100  CONTINUE 

BOUOO710 

0053 

GO  TO  800 

BOUOOT20 

0054 

810  1FIFQISI  .LT.  FQQI5II  GO  TO  800 

BOU00733 

0055 

A-AA 

BOU00740 

0056 

B-BB 

BOUOOT53 

0057 

C-CC 

BDU00760 

0058 

00  110  1-1,10 

BOUOOTTO 

0059 

FQIII-FQQIII 

80U00780 

0060 

110  CONTINUE 

80U00793 

0041 

800  CONTINUE 

BOUOOSOO 

0062 

MRITEI 6.6001  L2 

B0UJ0810 

0043 

400  FORMAT  1 LH1 , 25X, 'OPTIMAL  VALUE  FlMCTION  BOUNDS  HHEN  PARI", 12, 

6*1  IS  PERTURBED", /, 24X, 57I"P" 1, //I 

BOU00820 

B0U00830 

0064 

HR1TEI 6,6101  L2,E1,L2,FE1,L2,E2,L2,FE2 

80U00840 

0045 

410  FORMAT  12*," POINT  1  (UNPERTURBED  SOLUTIONI  t • ,/,2X, "PARI • , 12, 

6M-"  »E  15.T,10X,"F  1PAR1  *  ,12, *  il-*  ,E15.7,//,2R, "POINT  2  1  PERTURBED" 
8."  SOLUTIONI  I",/,2X,"PARI".I2,"II-",EI5.7,10X,"F(PARI",I2, 

»•»•", E15.7.///1 

80U00850 
B DU 00860 
80U008T0 
BOUOOSBO 

0046 

GO  TO  1910, 920,930,940, 950, 960I.NEXOP7 

80U00890 

c 

OUTPUT  HHEN  NEX0P7-1 

BDU00900 

0067 

910  HRITEI6.620I  S12.L2.812 

80U00910 

0048 

670  FORMAT  IPX, "LINE  PASSING  THROUGH  POINTS  1  ANO  2  ANO  OVER" 

B0U00920 

Figure  10. — Subroutine  B0UND. 


0069 

0070 

0071 

0072 

0073 

0074 


0073 

0076 


0077 

0078 

0079 

0080 

0081 

0082 

0083 

0084 

0083 

0086 

008  7 
0088 
0089 
0090 
0091 
0092 
0093 

0094 

0093 


0096 
039  7 

0098 

0099 

0100 

0101 

0102 

0103 

0104 

0103 

0106 

0107 

0108 

0109 

0110 

0111 

0112 

0113 

o:i4 

0113 

0116 


1371  NAT  I NG  F»  ',/,IX,S9C-'),//,llX.'F  .  *,115. 7,  •  •  PARC, 
$12,* 1  ♦  *,£15.7,///) 

MR  I  T  E  (  6, 6  30  I  SI.  1.2,81 

630  FORMATIZX.'LINE  UNDER  ESTIMATING  F*  AT  POINT  1  •  ./ 1 X  ,  361  I , 
i//.UX,'F  •  SEIS.T.'  •  PARI  • ,  12,  •  I  ♦  '.EIS.7.///I 
MR  I T  £  i 6,640 1  S2.L2, 82 

640  F0RMATI2X,'LINE  UNDER  E ST  I  MAT  INC  F*  AT  POINT  2* ,/. IX, 361 I. 
&//#llx,*F  ■  • , E 15.7, •  *  PARC,  12,  M  ♦  •  »E15,7,///I 
MR1TEI6.650I  A,L2,B,L2,C 

650  FORMAT I2X, *  QUADRATIC  ESTIMATION  OF  F»  THROUGH  POINTS  1  AND* 

*,■  2*,/,  IX, 501  l,  // ,1 IX,  • F  »  1 ,E1 5. 7, '  •  PARC. 12, 

Cl**2  ♦  •.E1S.7C  •  PARC,  12, •)  4  >,  E15.7  ,///! 

HR  1 1 E I  6,660 1  LZ 

660  FORMAT I2X,  • F*  BOUND  EVALUATION  AT  TEN  EOUIOISTANCE  POINTS', 

$•  BETMEEN  POINTS  1  AND  2  •  ,/,  IX,  7 1 C*'  I , //,  10X,  'PARC  ,  12,  •  C  , 
StOX, 'LOWER  BOUND', 10X,' UPPER  BOUND* ,1 OX, 'QUAD  ESTIMATE',/, 
»10X,TC-'I,9X,13C-'  ),8X,13C-'1.8X,I5C-'1,/I 
WRITE  16,6701 1 £( li,FL12l 1 1 .FU121 1 1 ,F0« 1 1 . I* 1 . 101 
MRITEI6,670IE2,FE2,FE2,FE2 
670  F0RMAT(10X,F7.3,7X,E15. 7,6X,E15.7,6X,E1S.7I 
RETURN 

C  OUTPUT  WHEN  NEXOP7«2 

920  MRITE(6,620IS12,LZ,812 
nR I T  E I  6,6601  LZ 

WRITE  I  6,921 1 (El 11.FU12I II ,  I >1,101 
WRITE! 6,921 1  E2.FE2 

921  FORMAT 110X,F7.3,28X,E15,7I 
RETURN 

C  OUTPUT  WHEN  NEX0P7-3 

930  WRIT£(6,«30I  SI,  LZ, 81 
WRIT  El  6,6601 S2,LZ*B2 
WRITE! 6,660 ILZ 

WRITE (6.931 1 1  El  1 1 .FL12I II ,I>1 , 101 
WRITEI6.931IE2.FE2 

931  F0RNAT(10X,F7.3,7X,E15.7I 
RETURN 

C  OUTPUT  WHEN  NEXOP7>4 

940  WR] TEI 6,941 I  S12.L2.B12 

941  FORMAT (2X,*L!NE  PASSING  THROUGH  POINTS  I  AND  2  AND  UNDER' 

ESTIMATING  F*  •,/,lX,59C-')»//,llXlr*F  -  ',E15.7.'  •  PARC, 
812, 'I  ♦  • ,E l 5. 7 ,///} 

WRITE  1  6,942 1  Sl.LZ.Bl 

942  FORMAT (2X, 'LINE  OVER  ESTIMATING  FP  AT  POINT  1  • , /IX, 361 I, 
8//  ,1 IX, 'F  -  ',£15. 7,'  *  PARC, 12, 'I  ♦  ',E15.7,///I 

WRITE! 6,9431  S2.LZ.B2 

943  F0«MATT2X,'L INF  OVER  ESTIMATING  F*  AT  POINT  2* , /, IX, 361 I , 
8//,UX,'F  -  '.£15.7,'  •  PARC, 12, •)  ♦ 

WRITE! 6,6501 A,LZ,B,LZ,C 
WR I T  E ( 6,660  I  L2 

WRITEI6, 670 l(E!tltFU12(ll,FL12lll»FQ(I|, 1*1,101 

WR  1 TE I  6,  670 1  E2.FE2.FE2.FE2 

RETURN 

C  OUTPUT  WHEN  NEXOP7-5 

950  WRITE! 6,942  I  Sl.LZ.Bl 
WRITE! 6,9431  S2.LZ.B2 
WRITE! 6,660  I  LZ 

WR1TE(6,921I(E1 I1.FL12! II,  1*1,101 

WRI TEI 6,921 1  E2.FE2 

RETURN 

C  OUTPUT  WHEN  NEX0P7*6 

960  WRITEI6,941IS12,LZ,812 
WRITE! 6,6601  LZ 

WRITE  16,931 1  I  El  I I.FU12! II,  1*1,101 

WRITE(6,931I  E2.FE2 

RETURN 

END 
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Figure  10. — continued 
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DIFF2.  Both  results  are  printed  out.  If  it  is  found  that  MATRIX 
creates  a  nonzero  element  below  the  main  diagonal  of  A  ,  then  a  switch 
is  set  to  cause  the  statement  ST0P  to  be  executed  in  MAIN.  Subroutine 
CHOKER  appears  as  Figure  11. 

6.4  C0NVRG 

After  each  iteration  of  the  algorithm  to  locate  the  minimum  of 
the  unconstrained  function  (the  W  function)  subroutine  C0NVRG  is  called 
to  determine  if  the  current  point  is  an  acceptable  approximation  of  a 
point  giving  the  minimum  value  of  the  W  function.  The  argument  N2  of 
this  routine  is  given  a  value  of  1  if  the  point  is  close  enough;  other¬ 
wise,  it  is  given  a  value  of  2.  C0NVRG  appears  as  Figure  12. 


6.5  DIFF1 


Subroutine  DIFF1  computes  the  first  derivatives  by  numerical 
differencing.  The  user-supplied  subroutine  RESTNT  is  called  2n  times 
for  each  gradient  evaluated  by  DIFF1.  For  the  function  f  ,  DIFF1  com¬ 
putes  the  1th  component  of  the  gradient  using  the  formula 


n  f(x°+0e  )  -  f(x°-0e  ) 

(Vf(x0)).  «  - ^ - ~ 


where  e^  is  a  vector  of  zeroes  with  a  1  in  the  ith  component  and 


0  is  a  small  number  whose  value  is  assigned  by  the  user, 
shown  as  Figure  13- 


DIFF1  is 


6.6  DIFF2 

Subroutine  DIFF2  computes  the  second  derivatives  by  numerical 
differencing.  The  matrix  of  second  partials  is  computed  using  central 
differencing.  The  differences  are  calculated  using  the  gradients  of  the 
function.  The  user-supplied  subroutine  GRAD1  is  called  2n  times  for 
each  matrix  of  second  partial  derivatives  evaluated  by  DIFF2.  DIFF2 
is  Figure  14. 
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DATE  -  6024? 


09/21/09 


SUBROUTINE  CHCKER 
C 

C  MARCH  1971 

C 

C  CHCKER  COMPUTES  AND  LIST  THE  FIRST  PARTIAL  DERIVATIVES  USING  GRAD1 
C  AND  THEN  USING  NUMERICAL  DIFFERENCING  (DIFF1I.  IF  REQUESTED  THE 
C  SECOND  PARTIAL  DERIVATIVES  ARE  COMPUTED  AND  LISTED  USING  MATRIX  AND 
C  DIFE2. 


CHE 00040 
CHFOOOSO 
CHE  00060 
CHE  00070 
CHE00040 
CHE  00090 
CHE 001 00 
CHE001 10 


0002 

IMPLICIT  REAL66C A-H.O-Z 1 

CHE00120 

0003 

REAL*4  XEP1.XEP2 

CHE 00130 

0004 

COMMON/ SHARE /XI 45 ltOELI4SI.AI45.45l ,N,N, NN.NPl .NM1 

CHE00140 

0005 

COMMON  /EQAL/  H,  HI,  M2 

CHE00150 

0006 

COMMON/E XPOPT/NEX OP l ,NEXQP2, NEX0P3, NEXOP4,NEXOP5,XEP1,XEP2 

CHE 00160 

0007 

NNZ-WN4NZ 

CHE 001 70 

0008 

DO  5  J-l.N 

CHE00180 

0009 

DELI Jl-l. 2345678 

CHE  00190 

0010 

5 

CONTINUE 

CHE  00200 

0011 

00  10  1-1, MMZ 

CHE 00210 

0012 

IN-I-1 

CHE 00220 

0013 

WRITE  16,1701  IN 

CHE 00230 

0014 

CALL  GRAD1  IINI 

CHE00240 

0015 

WRITE  16,1801  1  DEL! J! , J-l.NI 

CHE  00250 

0016 

CALL  DIFF1  IINI 

CHE00260 

0017 

WRITE  (6,1801  1  DELI J I, J-l ,NI 

CHE00270 

0018 

10 

CONTINUE 

CHE  00260 

C  •  •  • 

SOMETIMES  ONLY  FIRST  DERIVATIVES  ARE  TO  6E  CHECKEO 

CHE 00290 

0019 

IF  INEXOPl.LT. 41  GO  TO  160 

CHE00300 

0020 

WRITE  16,1901 

CHE00310 

0021 

00  150  I-l.MMZ 

CHE00320 

0022 

IN-I-l 

CHE00330 

0023 

WRITE  16,1701  IN 

CHE00340 

0024 

IT-2 

CHE00350 

0025 

DO  10  K-l.N 

CHE 00160 

0026 

00  20  J-l.N 

CHE00370 

0027 

20 

AIK.JI-O. 

CHF003B0 

0028 

30 

CONTINUE 

CHE 00390 

0029 

CALL  MATRIX  (IN.ITI 

CHE 00400 

00  30 

IF  1 IT.EQ.ll  GO  TO  150 

CHE 00410 

0031 

00  SO  K-2.N 

CHE00420 

0032 

KN1-K-1 

CHF00430 

0033 

00  40  J-I.KNl 

CHE 00440 

0034 

IF  IAIK.JI.E0.0.OI  GO  TO  40 

CHE 00450 

0035 

NEX0P1-5 

CHE 00460 

0016 

WRITE  16.2101  K, J 

CHE004T0 

0017 

GO  TO  60 

CHE  00460 

003B 

40 

CONTINUE 

CHE 00490 

0039 

50 

CONTINUE 

CHE00500 

0040 

60 

DO  90  K-l.N 

CHE 00510 

0041 

DO  70  J-K.N 

CHE 00520 

0042 

IF  IA|K, JI.NE.O.OI  GO  TO  60 

CHF00530 

0043 

TO 

CONTINUE 

CHE00540 

0044 

WRITE  (6,2201  K 

CHE  00550 

0015 

GO  TO  90 

CHE  00560 

0046 

60 

WRITE  16,2001  K,IAIK,JI,J-l,NI 

CHE  00570 

0047 

90 

CONTINUE 

CHE00580 

0048 

DO  110  K-1,N 

CHE  00590 

0049 

DO  100  2-1, N 

CHE  00600 

0050 

100 

AIK.JI-O. 

CHE 00610 

0051 

110 

CONTINUE 

CHE 00620 

0052 

WRITEI6.115I  IN 

CHE  006 10 

0053 

115 

FORMAT  IISHO  CALL  D1FF21 , 12, 1HI  1 

CHE  00640 

0054 

CALL  DIFF2  (INI 

CHF00650 

0055 

OO  140  K-l.N 

CHE00660 

0056 

DO  120  J-K.N 

CHE  00670 

0057 

IF  IAIK,JI.NE.OI  GO  TO  ISO 

CHE00680 

0056 

120 

CONTINUE 

CHE  00690 

0059 

WRITE  16,2201  K 

CHE 00700 

0060 

GO  TO  140 

CHE  00710 

0061 

110 

WRITE  (4,2701  K. (AIK, J». J-l.N! 

CHE  00720 

0062 

140 

CONTINUE 

CHE00T10 

0061 

150 

CONTINUE 

CHE  00740 

0064 

160 

CONTINUE 

CHE 00750 

0065 

RETURN 

CHE  00  760 

C 

CHE00770 

0066 

170 

FORMAT  1  20HQCMECKER . CONSTRAINT  NO.  ,  111 

CHE 00760 

0067 

160 

FORMAT  1 1 H0,24HCMECKER . 1ST  PARTIALS/I 1X,E20.8,E20.6,E20.6,E20 

.CHE  00 790 

18,E20.6,E20. 81 1 

CHFOOROO 

0066 

190 

FORMAT  <1M0,24HCMECKER . 2ND  PARTIALSI 

CHE  00610 

0069 

200 

FORMAT  I4HOROW,  IS  /I 1X,E20.6,c20.6,E?0.6,E20.6,E20.6,E20.6I 1 

CHE 00620 

0070 

210 

FORMAT  I1H  AI,I2,1H,»I2,10HI  .ME.  (.01 

CHE 00630 

0071 

220 

FORMAT  I4M  ROW.ll.UH  ALL  ZEROS.I 

CHE  00640 

0072 

END 

CHE 00650 

Figure  11. — Subroutine  CHCKER. 
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FORTRAN  IW  G  LEVEL  21  CONVRG  DATE  ■  B0262  09/21/37 


OOOl 

SUBROUTINE  CONVRG  IN1I 

CON00060 

c 

CONOOOSD 

c 

OCTOBER  1970 

CON 00060 

c 

CON00070 

c 

AFTER  EACH  ITERATION  OF  THE  ALGORITHM  TO  LOCATE  THE  MINIMUM  OF  THE 

CONOOORD 

c 

PENALTY  FUNCTION.  CONVRG  DETERMINES  IF  THE  CURRENT  POINT  IS  CLOSE 

C0N00090 

c 

ENOUGH  TO  THE  POINT  GIVING  THE  MINIMUM  VALUE  OF  THE  P  FUNCTION. 

cONOoino 

c 

N1  SET  EQUAL  TO  1  IF  MINtMUM  HAS  BEEN  FOUND. 

CONOOI 10 

c 

N1  SET  EQUAL  TO  2  IF  MINIMUM  HAS  NOT  BEEN  FOUND  ANO  TIME  IS  NOT 

UPCON00120 

c 

N1  SET  EQUAL  TO  )  OTHERWISE 

CONOOI 33 

c 

OOTT  SET  EQUAL  TO  (DEL  PI  ( INVERSE  (OELIOEL  PIIMOEL  PI  IN  OPT 

CON0016D 

0002 

IMPLICIT  REAL*8IA-H,0-2I 

CONOOISD 

0001 

REAL*6  RHOI N. RAT  10. EPS I tTHFTA0.XEPI.XEP2 

CONOOI 60 

0004 

CONMQN/SHARE/X 1651. DEL  1651, AI65, 651 .N.M.MN. NP 1  *  NMl 

CONOOI 70 

0005 

COMMON  /OPENS/  NTl,NT2.NT)rNT6.NT5,NT6,NT7,NT8.NT9,NT10 

CONOOIRO 

0006 

COMMON/ VALUE/F .G.PO.RSI GHA ,R J (SGI .RHO 

CONOOI9D 

0007 

COMMON/C RST/DFL X ( 651, DEL* 0(. , 5 1, RHOI N, RAT  10. EPS  I. THE TAO, 

C0N00200 

lRSlGl,Gl,Xll65l,X2l65l,X)i65l,XR2l65l,XRll65l,PRl, 

CON00210 

2PR2.Pl .F1.RJ1I90I.D0TT , PGR ADI  651. 01  AG (651, 

C0N0022D 

)  PREV5.ADELX.  NTCTR.  NUMINI.  NPHASE,  NSATIS 

CON00250 

0008 

COMMON/E XPOPT/NEXOPI,NEXOP2,NEXOP),NEXOP6,NEXOP5,XEP1,XEP2 

C0N0026D 

0009 

COMMON  /TSN/  NSWM 

C0N00250 

0010 

Nl-2 

CON00260 

OOtl 

IF  (MN.LE.tl  Ql-PO 

CON00270 

0012 

GO  TO  110, 20,501,  NT9 

C0N002B0 

001} 

10 

IFIDABSIDOrTI.LT. EPSII  GO  TO  TO 

CON00290 

0016 

GO  TO  60 

CON00500 

0015 

20 

lFIOABSIOOTTI.LT. (Pl-POI/S.UI  GO  TO  TO 

CONOO) 1 0 

0016 

GO  TO  60 

CON00520 

0017 

)0 

IF  (ADEL X.LT.EPS II  GO  TO  70 

CONOO 330 

0018 

60 

GO  TO  (50,601,  NSMM 

CQN00360 

0019 

50 

IF  (NN.LE.II  RETURN 

CQN0D350 

0020 

IF  (P0RXEP2  .LT.  Oil  GO  TO  75 

CQN00360 

0021 

WRITE  16,801 

CONOO 370 

0022 

GO  TO  70 

CONOO380 

002} 

60 

CALL  PUNCH 

CONOO 393 

0026 

WRITE  (6,901 

CONOOAOO 

0025 

Nl«) 

COND0610 

C 

CON 00620 

C 

FOUNO  THE  MINIMUM  TO  THE  SUBPROBLEN. 

CON00630 

0026 

RETURN 

C0N00660 

0027 

70 

Nl«l 

CON00650 

0028 

75 

01  •  PO 

CON00660 

0029 

RETURN 

CON 00670 

C 

CON00680 

00}0 

80 

FORMAT  UOOH  APPARENTLY  ROUNDOFF  ERRORS  PREVENT  A  MORE  ACCURATE 

DECON00690 

1TERMINAT ION  OF  THE  MINIMUM  OF  THIS  SUBPROBLEN. 1 

CON00500 

00)1 

90 

FORMAT  468H0****  TINE  IS  UP,  CALLING  EXIT  FROM  CONVRG.  ••••! 

CON 00510 

00)2 

ENO 

CON00520 

Figure  12. — Subroutine  CQNVRG. 
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0001 

SUBROUTINE  OtFFl  UNI 

D1F00040 

c 

OIF 00030 

c 

FEBUARV  1971 

DIF00060 

c 

OIF 00070 

c 

DIFF1  CONFUTES  THE  FIRST  DERIVATIVES  AT  NUNERICAL  DIFFERENCING. 

DIFOOUO 

c 

OIF00080 

c- 

-USER  CAN  CALL  FOR  DIFFERENCING  OF  SELECTED  FUNCTIONS 

DIF00100 

0002 

INFLICir  REAL*8( A-HtO-21 

DIFOOUO 

0003 

REAL*4  XEF1.XEF2 

DIF00120 

0004 

C0NNQN/SHARE/X14SI. DELI  45) >4143, 431, N.NtNN.NPl.NNl 

DIFOOUO 

0003 

C0NN0N/EXFQFT/NEX0Fl,NEX0F2,NEX0F3,NEX0F4tNEX0F5, XEF1, XEF2 

DIFOOUO 

0004 

C0NN0N/STIRX/XSTRt43l,XSSSI45l, DOLL  1451 

DIFOOUO 

OOOT 

00  10  J>1,N 

0IF00I60 

0008 

10 

XSTRI JI-XI  Jl 

OIF  001 70 

0008 

00  30  J»1.N 

DIFOOUO 

0010 

IF  IA.EQ.il  GO  TO  20 

DIFOOUO 

0011 

JN1* J— 1 

DIF00200 

0012 

XI JN1I-XSTRIJN1I 

DIF00210 

0013 

20 

X<JI-XS7R(J)»XEPl 

DIF00220 

0014 

CALL  RESTNT  (IN. 2221 

0 IF00230 

0013 

XI J)* XSTRI JI—XEP1 

0IF00240 

0016 

CALL  RESTNT  (IN.221I 

0IF00250 

ooir 

30 

DELI Jl>( 222-221 1/12. 4XEP1 I 

0IF00260 

0018 

XINI*XSTR(N) 

DIF00270 

0018 

RETURN 

DIF00280 

0020 

END 

DIF00290 

Figure  13. —Subroutine  DIFF1 
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Olf  F2 


DATE  ■  902A? 


OV/22/21 


0001 

SUBROUTINE  0IFF2  «IN! 

OIFOOOAO 

C 

OIFOOOAO 

c 

OCTOBER  1970 

OIFOOOAO 

c 

DIF00070 

c 

0IEF2  CONFUTES  THE  SECOND  DERIVATIVES  BY  NUNERICAL  DIFFERENCING. 

OIFOOOBO 

c 

OIF00090 

0002 

1  NFL  1C  IT  REALMI  A-H.O-Z! 

OIF  00100 

0003 

RE  AL  AA  XEPUXEP2 

0IF00110 

OOOA 

CONNON/SNARE/XIAD.OELI  A5l,A|AS,A5I.N,N,NN,NPl,NNl 

OIF00120 

ooos 

CONNON/EXPOFT/NEXQF1 «NE  X0P2. NExOPl. NE X0PA.NEX0P1.XEP1*  XEP2 

OIFOOIIO 

OOOA 

C0NN0N/ST1RX/XSTRIASI aXSSSlAlli DOLL  f All 

OIF  001  AO 

0007 

00  10  J*1,N 

OIFOOIIO 

0008 

10 

xsssui-xi.li 

OIF  00160 

0009 

oo  so  j*ijW 

01 F001 70 

0010 

IF  U.EO.ll  GO  TO  20 

DIF00180 

0011 

JN1-J-1 

OIF  001 90 

0012 

XUN1I-XSSSON1I 

OIF00200 

0013 

20 

XI JI-XSSSI JIaXEFI 

01 F  002 10 

001A 

CALL  GRA01  UNI 

0IF00220 

001 5 

00  10  1-1, N 

OIF 002 30 

001A 

10 

OOILIII-DELUI 

0IF002A0 

001? 

XI  J1«XSSSI  JI-XEP1 

0IF002S0 

0018 

CALL  GRA01  UNI 

0IF002A0 

0019 

00  AO  l*J,N 

OIF00270 

0020 

AO 

Al  J* I l-l DOLL  1 1 l-OELl 111/12. 9XEP1I 

0IF002B0 

0021 

10 

CONTINUE 

0 IF00290 

0022 

XINI-XSSSIttl 

0IF00103 

0021 

RETURN 

OIFOOIIO 

002A 

END 

0IF00120 

Figure  14. — Subroutine  DIFF£. 
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After  the  minimum  of  the  W  function  for  a  given  value  of  r  has 
been  located,  ESTIM  is  called.  When  the  minimum  of  the  W  function  for  a 
given  value  of  r  is  determined  a  subproblem  is  said  to  be  solved. 

ESTIM  performs  the  computations  to  estimate  the  Lagrange  multipliers  and 
make  the  first-  and  second-order  estimates  of  the  final  solution  of  the 
problem.  ESTIM  is  Figure  15. 

6 . 8  EVALU 

In  the  normal  phase,  EVALU  calls  the  user -supplied  routines  to 
evaluate  the  objective  function  and  the  constraint  functions  at  the 
current  point  x  .  As  the  constraint  functions  are  being  computed,  the 
code  performs  the  calculations  to  compute  the  penalty  terms  of  the 
penalty  function.  In  the  feasibility  phase  this  routine  puts  the  nega¬ 
tive  sum  of  the  violated  constraints,  which  is  the  objective  function  of 
the  entry  problem,  in  location  F  .  In  Location  F  the  current  value 
of  the  objective  function  is  stored  for  use  by  the  program.  After  F 
has  been  computed,  the  value  of  the  penalty  function  is  computed.  Also 
subroutine  EVALU  computes  a  value  for  the  dual  objective  function.  This 
value  is  only  correct  at  the  solution  to  a  subproblem  and  the  value  is 
stored  in  the  location  labeled  G  .  EVALU  appears  as  Figure  16. 

6.9  FEAS 

Subroutine  FEAS  determines  whether  the  starting  point  is  an 
interior  feasible  point  or  not.  If  the  variables  of  a  problem  are  sup¬ 
posed  to  be  nonnegative  and  some  of  the  components  of  the  starting  point 
vector  are  nonpositive,  then  those  components  are  changed  to  small  posi¬ 
tive  numbers.  A  modification  of  FEAS  must  be  made  if  it  is  desired  to 
change  this  number.  If  the  starting  point  does  not  satisfy  the  non¬ 
trivial  constraints  (the  constraints  other  than  the  nonnegativity  con¬ 
straints),  then  the  program  goes  into  the  feasibility  phase.  In  this 
phase  the  negative  of  the  sum  of  all  the  violated  inequality  constraints 
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FORTRAN  IV  C  LEVEL  21  ESTIN  DATE  •  80242  09/22/AS 


OOOl 

SUBROUTINE  ESTIN 

ESI 00040 

c 

E  ST00050 

c 

OCTOBER  1970 

f ST00060 

c 

EST00D70 

C  ESTIM  PERFORMS  THE  COMPUTATIONS  TO  ESTIMATE  THE  LAGRANGE  MULTIPLIERS 

EST00090 

C  ANO  MAKE  THE  FIRST-  AND  SECOND-ORDER  ESTIMATES  OF  THE  FINAL  SOLUTION 

EST00090 

C  OF 

THE  PROBLEM. 

EST00103 

0002 

IMPLICIT  REAL»8IA-H,0-2I 

E  ST  001 1 0 

0003 

RE AL *4  RHOIN.RATIOiEPSI tTHETAO 

EST00120 

00  OA 

COMMON/ SHARF /XI4SI. DELI  45 I.AI 45.451 ,N, M , MN, NP 1 . NM1 

EST00I30 

OOOS 

COMMON  /EOAL/  H,  HI.  M2 

E  ST  00140 

0006 

COMMON  /OPTNS/  NTl.NT2.NT3.NTA.NT5.NT6.NTTtNT8.NT9.NT10 

EST00150 

0007 

COMMON /VALUE/F iG.PO.RSIGMA.RJf  901 ,RHO 

EST00160 

0008 

COMMON/CRST/OEL  X(45I,DELX0(45I, RHOI N, RAT  1 0. EPS  I .THETAO, 

ESTOOITO 

lRSIGl,Gl,Xll45),X2(45),X3(4rl,XR2(45),XRII4SI.PRl, 

FST00190 

2PR2.P1.F1.RJ1I90) ,DOTT .PGR ADI A5I.DIAGIA51. 

EST00190 

3  PREV3.ADELX,  NTCTR,  NUMINI,  NPHASE.  NSATIS 

EST00200 

0004 

COMMON  /CONPAR/  NF1.NF2.NF3 

EST00210 

0010 

CALL  STORE 

EST00220 

0011 

210>RATIO**2 

EST00230 

0012 

29-RATIO 

E  ST  00240 

0013 

21-1. 0/29*1. 0/210 

E  ST  002  50 

001A 

22-21*1. /29PPJ 

EST00260 

001 5 

23-1./29PP3 

6  STOO?  70 

0016 

24-210*29 

EST00290 

0017 

25-2JPP3 

E  ST  00290 

0018 

26-1.0/(1210-1.01*129-1.011 

EST 00300 

0019 

27-1. /29 

EST00310 

0020 

28-1 . /I  29-1. t 

EST  00320 

0021 

RO-l.O/RHO 

E  ST  00330 

0022 

IF  ( NUMINI -21  150.80.10 

EST00340 

0023 

10 

WRITE  16,330) 

EST  00350 

002A 

PO-I PR 2-2A*PRl*25*P 11*26 

EST00360 

0025 

G-( R AT I0*G1-GR1 l/IRATIO-1.) 

EST00370 

0026 

DO  20  1-1,  N 

E  ST  00380 

0027 

20 

XIII-IXR 2(1) -2 A*XR1III*25*X 11111*26 

EST  00390 

0028 

NP-NPHASE 

EST00400 

0029 

NPHASE-A 

EST004I0 

0030 

CALL  EVALU 

E  ST  00470 

0031 

NPHASE-NP 

E  ST  00430 

00J2 

CALL  OUTPUT  121 

E  ST  00440 

C  CHECK  TO  SEE  IF  ESTIMATES  HAVE  CONVERGED 

E  ST00450 

0033 

GO  TO  170,30,70).  NPHASE 

EST00460 

003A 

10 

00  50  J-l.N 

EST  00470 

0035 

IF  IRJIJl)  40.50.50 

EST004B0 

0036 

40 

IF  (THETAOpRJIJ))  70.50.50 

EST  00490 

0017 

50 

CONTINUE 

E  ST  00500 

0038 

GO  TO  (70,70,60),  NTA 

EST0051 0 

0039 

60 

CALL  FINAL  INF3) 

EST00520 

OOAO 

TO 

CONTINUE 

E  ST  00530 

00  A 1 

SO 

NRITE  (6.3401 

EST00540 

00A2 

P0«(Z9*Pl-PRll*28 

EST00550 

00A1 

G-IR AT I0*G1-GR1)/ (RATIO-1.1 

EST00560 

OJA  A 

00  90  I-l.N 

EST00570 

00A5 

90 

XI 1 1-(  29PXU  D-XR1I1 11*28 

EST00580 

0046 

NP-NPHASE 

EST00590 

OOAT 

NPHASE-4 

EST00600 

00A8 

CALL  EVALU 

EST006IO 

0049 

NPHASE-NP 

EST  00620 

0050 

CALL  OUTPUT  12) 

EST00630 

C  CHECK  TO  SEE  IF  ESTIMATES  HAVE  CONVERGEO 

EST00640 

0051 

GO  TO  (140,100,140),  NPHASE 

EST  00650 

0052 

100 

DO  120  J-l.N 

EST00660 

005) 

IF  IRJIJl)  110,120,120 

E  ST  00670 

0054 

110 

IF  IRJIJI*THETAQI  140,120,120 

EST006B0 

0055 

120 

CONTINUE 

EST0069Q 

0056 

GO  TO  1140,130,1401,  NT 4 

EST00700 

0057 

130 

CALL  FINAL  (NF2I 

E  ST  00710 

0058 

1  AO 

CONTINUE 

EST00720 

0059 

150 

WRITE  16,350) 

E  ST  00730 

0060 

IF  IN)  180,180.160 

EST00740 

0061 

160 

DO  170  J-l.N 

EST00750 

0062 

170 

RJI JI-RHO/RJK  J) 

EST00760 

0063 

180 

IF  INZI  210,210,190 

EST00770 

0066 

190 

DO  200  J-l.NZ 

ESTOOTBO 

0065 

MNJ-M* J 

EST00790 

0066 

200 

RJI  MN  J  I-2.*RJ1(MNJ)*R0 

E  STOOBOO 

0067 

210 

60  TO  1220,240),  NT2 

ESTOOBIO 

0068 

220 

DO  230  l-l.N 

EST00A20 

0069 

2)0 

Kill -A  HO/ XI  III 

f STOOS30 

0070 

240 

CALL  OUTPUT  12) 

EST00840 

0071 

CALL  REJECT 

EST00B50 

0072 

IF  INUMINI-21  280,300,250 

EST00960 

007) 

250 

GO  TO  1280,310,260).  NTT 

E  ST  00970 

C 

SECOND  DROER  MOVE  FOR  NEXT  MINIMUM 

EST 00990 

00  74 

260 

on  270  I-l.N 

EST 00990 

0075 

270 

0ElXIII-2l*XIIII-22*XRlll l*2)*XR2ll 1 

EST00900 

Figure  15. — Subroutine  ESTIM 


T-434 


0076 

280 

PR2-PR1 

EST00910 

0077 

OR 2- GUI 

E  ST  00920 

0078 

PRl-Pl 

E ST 00930 

0077 

081-01 

EST00940 

0080 

00  290  I-l.N 

ESC 00950 

0081 

882(11-881111 

EST00960 

0082 

290 

881(11-81(11 

EST00970 

008) 

RETURN 

EST00980 

0084 

300 

CO  TO  I280«310*310>«  NTT 

E ST 00*90 

0085 

310 

DO  320  1-1. N 

estoiooq 

0086 

320 

DEL81 11-1811(1-8811 111*27 

EST010I0 

0087 

CO  TO  280 

f  ST01020 

C 

EST01039 

0088 

330 

PORN AT  1/26MO  2N0  ORDER  ESTIMATES  1 

EST01340 

0089 

340 

FORMAT  I/26H0  1ST  ORDER  ESTIMATES  1 

EST010S0 

0090 

350 

FORMAT  I/2SH0  LAGRANGE  MULTIPLIERS  1 

EST01Q60 

0091 

END 

EST01070 

Figure  15. — continued 
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FORTRAN  IV  6  LEVEL  21  EVALU  DATE  •  80242  04/28/17 


0001 

SUBROUTINE  EVALU 

EVA00043 

C 

EVA00050 

C 

OCTOBER  1970 

E VA00060 

C 

FVA00070 

C  IN 

THE  NORMAL  PHASE  EVALU  CALLS  THE  USER-SUPPLIED  ROUTINES  TO  E VALUATEEVA00080 

C  THE 

■  OBJECTIVE  FUNCTION  AND  THE  CONSTRAINT  FUNCTIONS  AT 

the  current 

EVA00090 

C  POINT.  IN  THE  FEASIBILITY  PHASE  THIS  ROUTINE  PUTS  THE 

NECATIVE  SUM 

OF  E VA00100 

C  THE 

!  VIOCATED  CONSTRAINTS  IN  LOCATION  F. 

EVA001 10 

0002 

IMPLICIT  REAL*8IA-H,0-ZI 

EVA00120 

0003 

REAL*4  RHOIN, RATIO. EPSI .THETAO 

EVA00130 

0004 

COMMON/SHARE /X| 45 1 • DELI 45 I.AI 45.451 .N.N.MN.NPl.NMl 

EVA00140 

0005 

COMMON  /EQAL/  H,  HI.  M2 

EVA00150 

0006 

COMMON  /QPTNS/  NTUNT2.NT3.NT4.NT5tNT6.NTT. NTS. NT9 

•  NTIO 

EVAO0I6O 

0007 

COMMON/VALUE/F.O.PO.RSISMA.RJISOI.RHO 

EVAOOITO 

0008 

COMMON/C RST/DEL XI 451 .DELXOf 45! ,RHO! N.RATIO.EPS 1 .THETAO. 

EVA00180 

1RSIG1.G1.X1I4S|,X2I45I.X3I45I.I:R2I45I.XR1|4SI.PR1. 

EVA0OI90 

2PR2.P1.F1.RJ1I90I.D0TT. PGR ADI 451 .01 AGI451 . 

EVA00200 

3  PREV3.ADELX.  NTCTR.  NUMINI.  NPHASE •  NSATIS 

EVA00210 

0009 

H-0.0 

EVA00220 

0010 

ASIGHA-0.0 

EVA00230 

0011 

F«0.0 

EVA00240 

0012 

NSATIS-2 

EVA00250 

0013 

GO  TO  110.100.190.2001.  NPHASF 

EVA00260 

C 

•1  FEASIBILITY 

EVA00270 

c 

•2  NORMAL 

EVA00280 

c 

•3  GUESS 

EVA00290 

c 

•4  ALL  FUNCTIONS  ARE  TO  BE  EVALUATED 

EVA00300 

C  FEASIBILITY 

EVA00310 

0014 

10 

GO  TO  120,401.  NT 2 

E  VA00320 

C  NON 

l-NE  GAT  IVIES  INCLUDED 

EVA00330 

00X5 

20 

DO  30  I'l.N 

EVA00340 

0016 

IF  IXIIII  260.260.30 

EVA00350 

0017 

30 

RSIGMA>RS!GMA-RHO*DLOGIXI  III 

E  VA00360 

0018 

40 

IF  IM.EQ.OI  GO  TO  90 

EVA00370 

0014 

00  80  J-l.M 

EVA00380 

0020 

CALL  RESTNT  IJ.RJIJII 

EVA00390 

0021 

IF  (RJlt JI.LE.O.OI  GO  TO  50 

EVA00400 

0022 

IF  IRJt Ji.GT.O.OI  GO  TO  60 

E VA00410 

C  VIOLATION  OF  A  PREVIOUSLY  SATISFIED  CONSTRAINT 

E VA00420 

0023 

GO  TO  260 

EVA00430 

0024 

50 

IF  IRJt JI.GT.O.OI  GO  TO  TO 

EVA00440 

C  ALL 

VIOLATEO  CONSTRAINTS  ADDED  INTO  OBJECTIVE  FUNCTION 

EVA00450 

0025 

F«F-RJIJI 

EVA00460 

0026 

GO  TO  80 

EVA00470 

0027 

60 

RS!GHA>RS!GNA-RHO*DLOGIRJI Jll 

EVA00480 

0028 

GO  TO  80 

EVA00490 

C  INDICATES  SATISFACTION  OF  CONSTRAINT! 10RN0REI 

EVA00500 

0029 

70 

NSATIS*! 

EVA00510 

0030 

RSIGMA*RS1GNA-RH0*DL0GIRJ( Jll 

E VA00520 

0031 

80 

CONTINUE 

EVA00530 

0032 

90 

CONTINUE 

EVA00540 

C  EQUALITIES  NOT  CONFUTED  IN  FEAS.  PHASE 

EVA00550 

0033 

PO«F»RSIGMA 

EVA00560 

0034 

G>F-RHO*DFLOAT|NI 

EVA005T0 

0035 

IFINT2.E0.il  G>G-RHO*DFLOATINI 

EVA00S80 

0036 

RETURN 

EVAD0590 

C  REGULAR  PHASE 

EVA00600 

0037 

100 

GO  TO  1110.1301,  NT 2 

EVA00610 

C  NON 

NEGATIVITIES  INCLUDED 

E VA00620 

0038 

110 

DO  120  1*1. N 

EVA00630 

0039 

IF  IXIIII  260,260.120 

EVA0Q640 

0040 

120 

RS1GMA*RSIGMA-RHO*OLOGIX|  III 

EVA 00650 

0041 

130 

IF  IM.EQ.OI  GO  TO  150 

EVA00660 

0042 

DO  140  J«1»N 

F VA006T0 

0043 

CALL  RESTNT  IJ.RJIJII 

EVA00680 

0044 

IF  IRJIJI.LE.O.OI  GO  TO  260 

EVA00690 

0045 

RSIGNA«RSIGNA-RH0*D10GIRJIJII 

E  VAOOTOO 

0046 

140 

CONTINUE 

E VA00710 

C  EVALUATE  AND  AOO  IN  EQUALITY  CONSTRAINTS 

EVA00720 

0047 

ISO 

CONTINUE 

EVA00730 

0048 

CALL  RESTNT  10, FI 

EVA00740 

0049 

IF  INI  1  180.180,160 

EVA00750 

0050 

160 

DO  170  I-1.NZ 

EVA00760 

0051 

J«I*N 

EVAOOTTO 

0052 

CALL  RESTNT  IJ.RJIJII 

EVAOOTSO 

C  ADO 

INTO  THIRO  TERN  OF  P  FUNCTION 

E  VA00T90 

0053 

H>H»IRJ{ JII9P2 

EVA00800 

0054 

1T0 

CONTINUE 

EVA00810 

0055 

H*H/RHO 

EVA00S20 

0056 

180 

PO*RSICNA»H 

EVA00810 

0057 

P0«F4P0 

EVA00840 

0058 

G*2. 6H-RH06DFL0AT IN 1 

E VA00850 

0059 

G*G«F 

EVA00860 

0060 

IFINT2.EQ.II  6>6-RH0*0FL0AT|NI 

EVA008T0 

C  DUAL  VALUE 

E VA00880 

0061 

RETURN 

E VA00890 

c  cucss  phase  not  coded 

EVA00900 

Figure  16.—  Subroutine  EVALU. 
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006? 

190 

RETURN 

EVA00910 

C - 

STRAIGHT  FUNCTION  EVALUATION  1  HAi N*FE AS  ONLY! 

EVA00920 

0063 

200 

CONTINUE 

EVA00910 

006* 

IF  IN.EO.OI  GO  TO  220 

EVA00960 

0065 

DO  210  1-1. H 

CVA00950 

0066 

CALL  RESTNT  II.RJIII) 

EVA00960 

006? 

210 

CONTINUE 

EVA009TO 

0068 

220 

CALL  RESTNT  10. FI 

E VA00960 

C  EQUAL  1  TV  CONSTRAINTS 

EVA 00990 

006? 

IF  INI  1  250.250.210 

EVA01000 

0070 

210 

DO  260  I-l.NI 

EVAOIOIO 

0071 

KZ-N*I 

E VA0I029 

00?2 

260 

CALL  RESTNT  IKZ.RUIK2II 

EVAOIOIO 

0073 

250 

RETURN 

EVA01060 

c  CONSTRAINTS  violated  not  so  before 

EVA01050 

0076 

260 

NSATIS-1 

EVA01060 

0075 

F0-10.E15 

EVAOIOIO 

0076 

RETURN 

EVAOIOIO 

0077 

ENO 

EVA01090 

Figure  16. — continued 
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is  labeled  for  use  as  the  objective  function  in  the  entry  problem.  Then 
the  routine  calls  B0DY  to  minimize  this  auxiliary  function  subject  to 
the  set  of  satisfied  constraints.  In  the  feasibility  phase  if  a  vio¬ 
lated  constraint  is  fortuitously  satisfied  an  immediate  return  to  FEAS 
is  made  and  a  new  entry  problem  is  begun,  including  the  newly  satisfied 
constraint  in  the  set  of  constraints  active  for  the  entry  problem.  Thus 
the  entry  problem  can  result  in  a  series  of  NLP  problems  being  partly 
solved. 


FEAS  also  contains  tests  that  indicate  when  a  problem  does  not 
contain  a  feasible  starting  point.  Such  information  is  printed  out,  and 
control  is  returned  to  MAIN  with  indicators  set  so  that  MAIN  begins  to 
try  the  next  NLP  problem  in  a  stack  of  problems.  When  the  entry  problem 
is  not  a  com  ex  programming  problem  the  test  indicates  only  that  the 
program  is  in  a  region  from  which  it  will  be  unable  to  locate  a  feasible 
starting  point.  The  user,  by  supplying  another  starting  point,  may 
cuase  the  algorithm  to  generate  another  sequence  of  points  that  does 
lead  to  a  feasible  starting  point.  Figure  17  shows  FEAS. 

6.10  FINAL 

Subroutine  FINAL  contains  the  test  used  to  determine  if  a  point 
satisfies  the  final  convergence  criterion  of  the  algorithm.  If  a  point 
does  satisfy  the  criterion  chosen  by  use  of  option  5  then  N2,  the  argu¬ 
ment  of  this  routine,  is  given  a  value  of  1;  otherwise  it  is  given  a 
value  of  2.  FINAL  is  called  following  the  computation  of  the  solution 
estimates,  which  are  made  after  the  solution  of  the  subproblem.  FINAL 
appears  as  Figure  18. 

6.11  GRAD 

The  gradient  of  the  W  function  is  given  by  the  formula 

m  m+p  2h  (x) 

VW(x,r)  *  Vf(x)  -  l  — y  .  Vg,  (x)  l  — - -  Vh.  (x)  . 

i-1  1  i*m+l  r  1 
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FEAS 


DATE  •  80242 


09/21/28 


0001  SUBROUTINE  FEAS  FEA00040 

C  FEAOOOSO 

C  AUGUST  I4T1  K A 00060 

C  FEAOOOTO 


C  FEAS  DETERMINES  WHETHER  THE  STARTING  POINT  IS  FEASABLE.  IF  IT  IS  NOT >  FEAOOOBO 
C  FEAS  LOOKS  FOR  A  FEASABLE  ONE.  IF  NONE  EXISTS.  A  MESSAGE  IS  PRINTEO  FEA00090 


C  AND  CONTROL  RETURNS  TO  NA1N.  FEA00100 


0002 

IMPLICIT  AEAL*BIA-H,0-2I 

FEAOOUO 

0001 

REAL *4  RHOINt RATIO.EPSI .THETAO 

FFA00120 

0004 

COMMON/ SHARE/X 1451 .OEL I4SI.AI4S.45I .N.N.NN.NPl.NMl 

FEAOOUO 

00  OS 

COMMON  /OPTNS/  NT1 .NT2.NT3.NT4.NTS.HT6. NT7.NTB.NT9.NT10 

FEA00140 

0006 

COMMON/ VALUE/F .G.PO.RSIGMA.RJf  901 .RHO 

FEA00150 

0007 

COMMON/C RST /OEL  XI 45 1 .DELX0I4SI. RHOINt RATIO.EPSI .THETAO. 
IRS IG1. Gl. XII 451, X21 4SI , X1I4SI, XR2I4SI.XR1 1451 .PR1, 
2PR2.P1.F1.RJ1 1901 fOOTr.PGRAOI45I.OIA?<45l. 

3  PREVS.ADELX,  NTCTR.  HUM INI,  NPHASE,  NSATIS 

FEAOOUO 

FEAOOtTO 

FEA001B0 

FEA00190 

0008 

NPHASE-1 

FEA00200 

0009 

GO  TO  <10,901.  NTZ 

FEA002I0 

0010 

10 

NF!X«1 

FEA00220 

0011 

DO  10  1-1, N 

FEA00230 

0012 

IF  (Xllll  20.20,30 

FEA00240 

0011 

20 

NFIX-2 

FEA00250 

0014 

XI 1 l-l.E-05 

FEA00260 

001S 

10 

CONTINUE 

FEA00270 

0016 

GO  TO  <50,401,  NF IX 

FEA002B0 

0017 

40 

NPHASE -4 

FEA00290 

0018 

CALL  EVALU 

FEA00300 

C  JUST  GET  ALL  CONSTRAINTS  EVALUATEO 

FEAOOUO 

0019 

NPHASE-1 

FEA00320 

0020 

WRITE  <6,1301 

FEA00330 

0021 

CALL  OUTPUT  <21 

FEA00140 

0022 

50 

IF  INI  90,90,60 

FEA00350 

0023 

60 

00  70  1*1, M 

FEA00360 

0024 

IF  IRJIIII  100,100.70 

FEA00370 

002  5 

70 

CONTINUE 

FEA00180 

0026 

80 

CALL  TINEC 

FEA00390 

0027 

WRITE  <6,1401 

FEA00400 

0028 

G-0.0 

FEA004I0 

0029 

CALL  RESTNT  IO.FI 

FEA00420 

0030 

CALL  OUTPUT  121 

FEA00430 

0031 

90 

RETURN 

FEA00440 

0012 

100 

CALL  BOOV 

FEA00450 

0033 

IF  1 NPHASE  .EO.  51  RETURN 

FEA00460 

0034 

00  110  l-l.N 

FEA00470 

003S 

IF  IRJIIII  120,120,110 

FEA004BO 

0036 

110 

CONTINUE 

FEA00490 

0017 

GO  TO  80 

FEA00500 

0038 

120 

WRITE  <6.1501 

FEA00510 

C  TO 

INDICATE  TO  MAIN  TO  START  ON  NEXT  PROBLEM. 

FEA00520 

0019 

NPMASE-5 

FEA00530 

0040 

C 

GO  TO  90 

FEA00540 

FEA00550 

0041 

130 

FORMAT  < 1H0, 2X, 4SHMADE  VIOLATED  NON-NEGATIVITIES  SLIGHTLY 
11 

POSITI VEFEA00S60 
FEA00570 

0042 

140 

FORMAT  ISIH066P6PTHE  FEASIBLE  STARTING  POINT  TO  BE  USEO  IS  ...1  FEA00580 

0043 

ISO 

FORMAT  I1X.89HTHIS  PROBLEM  POSSESSES  NO  FEASIBLE  STARTING 
I  ILL  LOOK  FOR  DATA  FOR  NEXT  PROBLEM.  ) 

POINT,  WFEA00590 
FEA00600 

0044 

END 

FEA00610 

Figure  17. — Subroutine  FEAS 
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FINAL 


OATt  •  B02A2 


SVII/U 


FINOOOAQ 


0001 

SUBROUTINE  FINAL  IN2I 

FINOOOSO 

c 

F INOOOA3 

c 

OCTOBER  1970 

F IN00070 

c 

FINOOOSO 

C  FINAL  CONTAINS  THC  TESTS  USED  TO  OETERNINE  WHETHER  A  FOINT  SATISFIES 

FIN00090 

C  THE  FINAL  CONVERGENCE  CRITERION  CHOOSER  TO  OETERNINE  IF  THE  NLF 

F INOOIOO 

C  FR08LEM  HAS  REEN  SOLVED. 

FIN00I10 

C  N2 

SET  EQUAL  TO  1  IF  CONVERGENCE  CRITERION  IS  SATISFIED. 

F IN00120 

C  N2 

SET  EQUAL  TO  2  OTHERWISE • 

FINOOUO 

0002 

INPLICIT  REAL*StA-H,0-ZI 

FIN001A0 

0001 

REALAA  RHOIN. RATIO, EPSI.THETAO 

FIN001S0 

00  OA 

CONNON/SHARE/XIASI,DEltAS>,AIAS,ASl,N,N,NN,NPI,NNl 

FINOOIAO 

0001 

COHNON/VALUE/F.G, POfRSIGHA, RJI90I ,AHO 

FIN00I70 

0006 

CONNON/CRST/OELXIASI.OELXOI AS), RHOIN, RATIO, EPSI.THETAO. 

FINOOIAO 

IRSIGl.GltXll ASI,X2IASI,X1I*SI,FR2(ASI,XR1IASI,PR1, 

FINOOIAO 

2PR2,P1,F1,RJ1I90>,D0TT,PGRAD|A5>,DIAGIAS», 

FIN00200 

1  PREVl.AOELX,  NTCTR,  NUHINI,  N*HA$E,  NSATIS 

F IN00210 

0007 

GO  TO  110,20,101,  NTS 

FIN00220 

0008 

10 

EPS!L*OABStF/G-l.l 

FIN00210 

0009 

IF  IEPSIL-THETAOI  50, SO, TO 

F IN002A0 

0010 

20 

IFIOABSIRSIGNAI-THETAOI  SO, SO, TO 

F INO02S0 

0011 

10 

IF  INUWINI-ll  SO, AO, AO 

F IN002A0 

0012 

AO 

PESr>PRl-(PRl-POI/li.-l./SQAT(RAUOII 

FIN00270 

0011 

EPSIL«0ABSIPEST/G-1.I 

F1N002BO 

001 A 

IF  IEPSIL-THETAOI  30,70,70 

FIN00290 

0011 

SO 

N2«l 

FIN00100 

001A 

GO  TO  ISO.AOI,  NT A 

FINOOUO 

0017 

AO 

CALL  PUNCH 

FIN00120 

ooia 

GO  TO  BO 

FINOOUO 

0019 

70 

N2*2 

FINOOIAO 

0020 

80 

RETURN 

finooiso 

0021 

END 

FINOOIAO 

Figure  18. — Subroutine  FINAL 
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Subroutine  GRAD  calls  the  user-supplied  subroutine  GRAD1  to  compute 
Vf(x)  ,  Vg^(x)  ,  and  Vh^(x)  and  performs  the  computations  to  evaluate 

VW(x,r)  .  The  negative  of  the  gradient  of  the  U  function  (i.e.,  -Vw) 
j  left  in  the  array  DELXO  when  GRAD  returns  control  to  the  calling  rou¬ 
tine.  When  the  argument  of  GRAD  has  a  value  of  2  this  is  all  that  is 
done.  However,  when  it  has  a  value  of  1,  GRAD  also  does  part  of  the 
computations  needed  to  evaluate  the  matrix  of  second  partial  derivatives 
of  W  at  x  .  Subroutine  GRAD  is  shown  in  Figure  19. 

6.12  INVERS 

When  Newton's  method  is  used  to  minimize  the  W  function  for  a 
given  value  of  r  ,  it  maps  the  negative  of  the  gradient  of  W  with  the 
inverse  of  the  matrix  of  second  partial  derivatives  of  W  evaluated  at 
x  .  A  positive  definite  matrix  will  always  give  a  vector  along  which 
the  value  of  W  will  initially  decrease.  That  is,  for  iteration  i  ,  a 

move  is  made  along  the  vector  S1  ,  given  by  the  formula 

S1  =  -[V2W(x,r)3_1  VW(x,r)  . 

This  is  equivalent  to  solving  the  set  of  simultaneous  linear  equations 

[V2W(x,r)]Si  =  -VW(x,r) 

i  *  2 

for  S  ,  where  V  W  and  VW  are  respectively  the  matrix  of  second 
derivatives  and  gradient  vector  of  W  evaluated  at  (x,r)  . 

Subroutine  INVERS  solves  this  set  of  equations  for  S1  using  an 
L-U  decomposition  method  (the  Crout  procedure) .  If  it  is  determined 
that  the  matrix  of  second  partials  is  not  positive  definite,  a  different 

procedure  is  used  to  obtain  a  direction  S*  .  A  complete  discussion  of 
this  procedure  is  given  on  page  167  of  Fiacco  and  McCormick  [10].  When 


*ln  the  program  the  vector 


S1  is  called  DELX. 
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0001 

SUBROUTINE  GRAD  MSI 

CRA0D0*0 

c 

graoooso 

c 

OCTOBER  1970 

GRA00060 

c 

GRA00070 

C  GRAO  COMPUTES  THE  GRADIENT  OF  THE  PENALTY  FUNCTION  A  HO  THE  OUTER 

GRAOOORO 

C  PRODUCT  FACTORS  OF  THE  MATRIX  OF  SECONO  PARTIAL}  OF  P. 

GRA00090 

C  IF  US«1I  ACCOM.  MATRIX  OF  2ND  PART  1 ALS  IFIIS*2I  OONT 

CRA00100 

0002 

IMPLICIT  REAL»8I A-H.0-21 

GRA001 10 

0003 

REAL**  RHOIN.RATIOtEPSliTHETAO 

GRA00123 

000* 

C0MM0N/SHARE/XI*5 It  DEL <*SI.AI *$,*$» ,N,N,MN,NP1,NN1 

GRAOOI 30 

000$ 

COMMON  /EOAL/  H.  HI.  M2 

GRA001*0 

0006 

COMMON  /OPTNS/  NT1 ,NT2«  NT3,NT*«NT5, NT6.NT7 .NTS .NT9. NT10 

GRAOOI SO 

0007 

COMMON/VALOE/F .G.P0.RSIGMA.RJI90I.RH0 

GRA00160 

OOOS 

COMHON/CRST/OELXI*$l . OELXO! *5 I.RriOINt RATIO. EPS 1 .THETAO. 

GRAOOI TO 

lRSlGl.Gl,Xlt*S),X2l*5),X3t*'l .XR2I*SI .XR1IA5I ,PRl, 

GRA30180 

2PR2.P1 .F1.RJ1I90I .OOTT . PGRAOIASI . 01 AGI65I t 

GRAOOI 90 

3  PREV3.A0ELX,  NTCTR.  NUMINI.  NPHASE.  NSATIS 

GRA00200 

0009 

GO  TO  110.301.  IS 

GRA00210 

0010 

10 

00  20  1*1. N 

GRA00220 

0011 

00  20  J*1,I 

GRA00230 

0012 

20 

Atl.JI-O. 

GRA002*0 

0013 

30 

DO  *0  I-l.N 

GRA002S0 

001* 

*0 

DELXOIH-0. 

GRA00260 

C  THIS  SECTION  WORKS  CORRECTLY  IN  FEASIBlllTV  PHASE  AS  HELL  AS  NORMAL 

PHGRA00270 

001$ 

GO  TO  ISO. 801.  NT2 

GRA00280 

0016 

SO 

00  70  1*1 »N 

GRA00290 

0017 

0ELX0II1— RHO/XIII 

GRA00300 

0018 

GO  TO  160.701,  IS 

GRA00310 

0019 

60 

AII,ll*l-DELXOIII/XIIII 

GRA00320 

0020 

70 

CONTINUE 

GRA00330 

0021 

80 

CONTINUE 

GRA003*0 

0022 

IF  IN.LE.OI  GO  TO  180 

CRA00350 

0023 

00  170  K*l,N 

GRA00360 

002* 

CALL  GRA01  IK) 

GRA00370 

002$ 

IF  IRJIKI.GT.O.OI  GO  TO  110 

GRA00380 

C  ALL 

VIOLATED  CONSTRAINT  GRADS  ADDED  TO  OBJ.  FUNCTION 

GRA00390 

0026 

00  100  l«l,N 

CRAOO*0O 

0027 

IF  IOELIIII  90.100.90 

GRA00A10 

0928 

90 

OELXOI I l-DELXOI I I-DELl I 1 

OR AO 0*20 

0029 

100 

CONTINUE 

G* A 00* 30 

0030 

GO  TO  170 

GRA00660 

0031 

110 

TT-RHO/RJIK) 

graooaso 

0032 

00  160  1*1, N 

CRA00A60 

0033 

IF  IOELIIII  120,160,120 

GRA00670 

C  IF 

DELI  1 1-0  SKIP  ALL  THE  FOLLOWING  COMPUTATION  INVOLVING  •  BV  DELHI 

CRA00*80 

003* 

120 

T-TT6DELI11 

GRA00A90 

003$ 

OELXOI II*0ELX0III-T 

CRAOOSOO 

0036 

GO  TO  1190,1601,  IS 

CRA00510 

0037 

130 

T-T/RJIKI 

GRA00S20 

0038 

00  ISO  JJ*1,I 

GRA00S30 

0039 

IF  IOELIJJII  1*0,150,1*0 

GRA00$«0 

00*0 

1*0 

AII,JJI«A{|,JJI*T*OELIJJ) 

GRAOOSSO 

00*1 

ISO 

CONTINUE 

GRA00560 

00*2 

160 

CONT INUE 

GRA00570 

00*3 

170 

CONTINUE 

GRAOOSSO 

C  EOUALITY  CHANGES  FOR  GRAD 

GAAOOS90 

00** 

180 

IF  INZ.LE.OI  GO  TO  250 

GR A 00600 

00*$ 

GO  TO  I250,190«250l,  NPHASE 

GRA00610 

00*6 

190 

RO-2./RHO 

GRA00620 

00*7 

DO  2*0  J*1,NZ 

GRA00630 

00*8 

K*M*J 

GR*006*0 

00*9 

CALL  GRAOl  IKI 

GRAOOSSO 

OOSO 

TT*RO*RJ(KI 

GRA00660 

00$  1 

00  230  I-ItN 

GRA00670 

U0$2 

IF  IOELI I I.EQ.O.OI  GO  TO  290 

GRA00680 

OOS3 

OELXOI 1 l*DELXOI 1 )*DELI 1 l*TT 

GRA00690 

OOS* 

GO  TO  1200,2301,  IS 

GRA00700 

OOSS 

200 

T*R3*0EL 1 I 1 

GRA00710 

00  $6 

DO  220  JJ«I,I 

GRA00720 

OOS  7 

IF  IOELIJJII  210,220,210 

GRA00730 

0038 

210 

Al 1, JJI-AI I, JJI*T*OEL< JJI 

GRAO0T*O 

00S9 

220 

CONTINUE 

GRA00750 

0060 

230 

CONTINUE 

CRA00760 

F061 

2*0 

CONTINUE 

GR ROOT 70 

0062 

250 

GO  TO  1260,2801,  IS 

GRAOOTBO 

0063 

260 

DO  270  1*1, N 

GR  A 00 790 

006* 

270 

01 AGI I l*A( 1,11 

GRA00800 

006$ 

280 

GO  TO  1290,330,2901,  NPHASE 

GRA00810 

C  LEAVES  NEGATIVE  SRAOIENT  IN  OELP 

GRA00820 

0066 

290 

DO  300  I»l, N 

GRA00830 

0067 

300 

OELXOI  II —OELXOI  II 

GRAOOSSO 

0068 

310 

AOELXO. 

GRAOOSSO 

0069 

DO  320  1*1, N 

GRA00660 

0070 

320 

ADEL  X-AOEL X, OELXOI 1 1 6*2 

GRA00870 

0071 

AOELX-DSORTIAOELXI 

GRA00880 

0072 

RETURN 

GRA00890 

0073 

330 

CALL  GRAOl  101 

GRA09900 

007* 

DO  3*0  1*1, N 

6RA009I0 

007$ 

9*0 

OELXOI 1 l—DELXO 1 1 l-OEL 1 1 1 

GRA00920 

C  LEAVES  THE  NEC.  GR AO  OF  P  IN  OELXO 

GRAOOOSO 

0076 

GO  TO  310 

6RA00940 

007  7 

END 

GRAOOOSO 

Figure  19. — Subroutine  GRAD. 
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this  occurs  the  program  prints  out  "ORTHOGONAL  MOVE."  This  also  warns 
the  user  that  the  problem  is  probably  not  a  convex  program.  Subroutine 
INVERS  appears  as  Figure  20. 

6.13  LMULT 

Subroutine  LMULT,  called  by  subroutine  SENS,  calculates  the 
Lagrange  multiplier  sensitivities  according  to  steps  5  and  6  of  Algo¬ 
rithm  2.1  (Section  2).  LMULT  is  given  in  Figure  21. 

6.14  MAIN 

MAIN  is  the  program  that  initiates  the  SENSUMT  algorithm  to  solve 
a  nonlinear  programming  (NLP)  problem;  it  is  not  a  subroutine.  The  in¬ 
put  of  parameters  and  options,  as  well  as  a  starting  point,  is  done  in 
MAIN,  and  the  call  to  READIN  (a  user-supplied  subroutine)  is  made  to 
allow  the  user  to  read  in  data  needed  to  evaluate  his  objective  function 
and  constraint  functions.  Starting  point  data  or  blank  cards  should 
always  be  supplied  by  the  user  because  starting  point  data  cards  are 
always  read  in  MAIN.  Subroutine  FEAS  is  called  to  obtain  a  feasible 
starting  point  by  making  use  of  the  SENSUMT  algorithm  to  solve  the  entry 
problem  if  the  user-supplied  point  is  not  feasible.  The  actual  solution 
of  the  NLP  problem  using  the  SENSUMT  algorithm  is  supervised  by  subrou¬ 
tine  B0DY.  The  calls  to  SENS  (sensitivity  subroutine)  and  B0UND  (bound 
subroutine)  are  also  done  in  MAIN.  When  calculating  bounds,  after  the 
solution  and  sensitivity  analysis  of  the  unperturbed  problem,  MAIN  re¬ 
iterates  SENSUMT  to  solve  the  subject  problem  with  the  perturbed  data. 
Perturbations  and  readjustments  in  the  problem  data  are  done  in  PERT 
(perturbation  subroutine),  on  call  by  MAIN.  Subroutine  MAIN  is  shown 
as  Figure  22. 

6.15  0PT 

The  purpose  of  subroutine  0PT  is  to  obtain  a  §  >  0  such  that 
W  (x  +  OS)  is  a  minimum  with  respect  to  9  along  the  vector  x  +  0S  , 
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DATE  •  >0242  12/20/20 


0001  SUBROUTINE  INVERS  (NSHEI  INV00060 

C  INV00050 

C  OCTOBER  19  70  INV00060 

C  INVOOOTO 

C  INVERS  SOLVES  THE  SET  OF  EQUATIONS  FOR  THE  MOVE-VECTOR  USING  THE  INVOOORO 

C  CROUT  PROCEDURE.  IF  THE  MATRIX  IS  NUT  POSITIVE  0EF1N1TC.  A  OlFfeUBN T  INV0OO9O 
C  METHOO  IS  USED.  INV00103 

C  INVOOliO 

(••••PERFORMING  A  L-U  DECOMPOSITION  OF  THE  MATRIX  A.  TAKING  AOVANAGE  OF  INV00123 
(••••THE  SYMMETRY  OF  THE  A  MATRIX.  INV00130 

(••••IF  A  NON-POSITIVE  PIVOT  CANDIDATE  IS  GENERATED.  THEN  MCCORMICK, S  INV00160 

C* •••PROCEDURE  IS  USEDISEE  PP.  167-168  IN  FIACCO  AND  MCCORMICKI.  INV001S0 

C  I NV 001 AO 

(••••IF  NSME  >1  NORKING  HITH  A  NEM  A  MATRIX.  IF  NSME*  2  USING  PREVIOUS  INVOOITO 
(••••A  MATRIX,  BUT  HAVE  A  NEN  RIGHT-HAND-SIDE.  INVOOIBO 

C  1NV00190 

(••••MINV  IS  THE  NUMBER  OF  NON-POSITIVE  PIVOT  (ANDI DATES  GENERATED.  INV00200 


0002 

IMPLICIT  RE AL*8 I A-H, 0-2 1 

INV00210 

0003 

REAL*6  RH01N, RATIO, EPSI .THETAO, XEPl ,XEP2 

1 NV00220 

0006 

COMMON/SHARE/X ! 65) , DEL! 65 !#At 65,65) ,N, M,NN, NP1 ,NM1 

INV00230 

ooos 

COMMON  /OPT NS/  NT1,NT2,NT3,NT6,NT5,NT6,NT7,NT8,NT9,NT10 

I NV00260 

0006 

COMMQN/CRST/DELXI 65), OELXOI 65), RHOIN, RATIO, EPSI, THETAO, 

INV00250 

lR$IGl,Gl,Xl|65l,X2I65l,Xll65l,XR2l6SI,XRll65),PRt, 

INV00260 

2PR2 ,P1,F1,RJ1I90I ,DOTT , PGRADlASIaDI AGI65I , 

INV00270 

3  PREV3.A0ELX,  NTCTR,  NUNIN1.  NPHASE,  NSATIS 

INV00280 

0007 

COMMON/EXPOPT/NEXOP1,NEXOP2,NEXOP3,NEXOP6,NEXOP3,XEP1,XFP2 

INV00290 

0008  ' 

DIMENSION  8165) 

I NV00300 

0009 

EQUIVALENCE  IB.OELXI 

INV00310 

0010 

GO  T0120*  170) ,  NSME 

I NV00320 

0011 

20 

NINV*0 

INV0033O 

0012 

IF  IAIl.il)  60.30,50 

INV00360 

0013 

30 

N1NV-1 

INV00350 

0016 

GO  TO  70 

INV0036Q 

0013 

60 

NINV-l 

1NV00370 

0016 

50 

All. 11*1. /All, 1) 

INV00380 

0017 

DO  60  1-2, N 

I NV 00390 

0018 

60 

All, I) -A 11, I 1*111,1) 

INV00600 

0019 

70 

00  160  4-2, N 

INV00610 

0020 

JMl-J-1 

INV00620 

0021 

7-0. 

INV00630 

0022 

DO  90  I-1.4N1 

INV00660 

0023 

IF  IAI 1, J) |  80,90.80 

INV00650 

0026 

80 

T-T*AIJ,II*AII,J) 

I NV00660 

0023 

90 

CONTINUE 

INV00670 

0026 

Al J.JI-AI J.JI-T 

I NV 006 80 

0027 

IF  lAIU.JI)  110,100,120 

INV 00690 

0028 

100 

NINV-NINV»1 

INV00500 

0029 

GO  TO  170 

■ NV00510 

0030 

HO 

NINV-NINVH 

INV00520 

0031 

120 

Al J, JI-1./A1 J.J) 

INV 00330 

0032 

IF  IU.EQ.N)  GO  TO  170 

INV00560 

0033 

JP1-J*1 

INV 00550 

0036 

DO  130  L-JPl.N 

INV00560 

0035 

T»0. 

INV005T0 

0036 

DO  160  1-1, JMl 

INV 00580 

0037 

IF  IAI1.JI)  130.160.130 

INV 00 590 

0038 

130 

T-T»AIL,))MII,J> 

INV00600 

0039 

160 

continue 

1NV00610 

0060 

All,  JI-All,J»-T 

INV 00620 

0061 

A 1 J.L l-AIL,4l*AIU, J) 

INV00630 

0062 

ISO 

CONTINUE 

1 NV00660 

0063 

160 

CONTINUE 

I NV00650 

0066 

170 

CONTINUE 

I  NV  00  A  60 

0063 

IF  ININVI  180,180,290 

INV00670 

0066 

180 

Blll*Blll*All,tl 

I NV00680 

0067 

00  210  J-2.N 

INV 00690 

0068 

T-0. 

INVOOTOO 

0069 

JN1-4-1 

INV00710 

0030 

DO  200  I-1.4NI 

INV 00 720 

0031 

IF  1 A| J, || |  190,200,190 

INV 00730 

0032 

190 

T-T«AIJ, 11*8111 

INV00760 

0033 

200 

CONTINUE 

INV0075 0 

0036 

81 J 1*1 81 J l-T l*AI J, J I 

1 NV00760 

0033 

210 

CONTINUE 

INV00770 

0036 

DO  260  I-1,NN1 

INVOOTBO 

0037 

NMK-N-I 

INV00790 

0038 

1 

DO  230  J-l, 1 

I NV00800 

0039 

L-NP1-J 

iNvooaio 

0060 

IF  lAlNMK.LII  220,230,220 

INV 00820 

0061 

220 

BINM«)-BINNKI-AINNK,L)*BILI 

INV 00830 

0062 

230 

CONTINUE 

INV00860 

0063 

260 

CONTINUE 

INV00830 

0066 

230 

CO  TO  1280,2601,  NTS 

INV00860 

0063 

260 

URITE  16,6301 

INV00870 

0066 

MRITE  16,6201  IDELXOII > ,1-1 ,NI 

INV 00880 

0067 

270 

NRITE  16,6601 

INV 00890 

0068 

MRITE  16.6201  IDELXI 1 1 , i-1 .N) 

1 NV00900 

Figure  20.— Subroutine  INVERS 


006* 

260 

HE  TURN 

1 NV00810 

t - 

COMPUTE  ORTHOGONAL  MOVE 

INV00820 

0070 

260 

CUNT INUE 

1 NV00930 

0071 

00  >50  1 1*1 »  N 

INV00640 

0072 

(•N-II*1 

1 NV00950 

0073 

IF  (AM. Ill  >10. >00, 320 

INV00960 

007* 

>00 

Bl  11*0.0 

INV 00670 

0075 

GO  TO  >50 

INV006B0 

0076 

>10 

6(11*1.0 

INV 00980 

0077 

GO  TO  »0 

INVOIOOO 

007B 

>20 

•(11*0.0 

INVU1010 

007* 

»0 

IPI*1»1 

INV01020 

ooao 

IF  (IP1.GT.NI  GO  TO  >50 

1 NV01030 

0081 

DO  340  J"IP1,N 

INV01040 

0082 

>40 

B(I>*B(II-AII,JI*B(JI 

I NV01050 

0083 

>50 

CONTINUE 

1  NV01060 

0084 

GO  TO  >60 

INV01070 

C—  CHECK  MAYBE  00  OIFF  FOR  P.S.D. 

1 NV01060 

0086 

360 

2C2-0.0 

1  NV01090 

0086 

00  370  1*1, N 

INV01100 

0087 

370 

ZC2*ZC2*DELX0I  I !*B(  1 1 

INVOUIO 

0088 

IF  IZC2I  >80,400.400 

INV01I20 

0084 

380 

00  380  1*1. N 

INV01D0 

0090 

390 

BUI— Bill 

INV01140 

00*1 

400 

WRITE  16,4501 

INV01150 

l  MCC 

ZANGMIU.  ONE  MOO 

INV01160 

0092 

IF  IMEXOP2.NE.2I  GO  TO  250 

INV0U70 

00*3 

DO  410  K*1,N 

INV0U80 

00*4 

410 

S(KI*B(K)«0ELXO(K) 

INV01190 

00*5 

GO  TO  250 

1NV01200 

c 

INV01210 

0096 

420 

FORMAT  17E1T.BI 

1NV01220 

00*7 

430 

FORMAT  ( 1H0.6X, 12H0EL  P  VECTOR! 

INV01230 

0098 

440 

FOANATOHO, 6X ,24HS ECONO  ORDER  MOVE  VECTORI 

INV0I240 

0099 

450 

FORMAT  ( 1H0.&X,  15  HORTHOGONAL  MOVE! 

INV01250 

0100 

END 

INV01260 

Figure  20. — continued 
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FORTRAN 

IV  G  LEVEL 

21 

LMULT 

DATE  *  802A2 

12/08/23 

LNU00040 

0001 

SUSK0JT1NC 

LMULTI  INO.OELNU.PEM.OUI 

t  *1000050 

C 

LNU00060 

c 

1  MARCH  1976 

LHU00070 

c 

LNU00093 

C  SUBROUTINE  LMULT  IS  USED  TO  ESTIMATE  THE  PARTIAL  DERIVATIVES  OF  THE  LMU00090 
C  LAGRANGE  MULTIPLIERS.  IT  IS  CALLED  BV  SENS  WHEN  NEX0P6  •  1.  USING  LNUOOIOO 
C  THE  DIFFERENTIABILITY  OF  Xlllt  IT  TAKES  ADVANTAGE  OF  THE  RELATIONS  LNUOOIIO 
C  Ulll  >  RHO/GI 1 1**2  AND  tfUl  »  2»H(  JS/RHO,  DIFFERENTIATING  THEM  WITH  LNU00120 
C  RESPECT  TO  THE  PARAMETERS.  THIS  SUBROUTINE  MAS  CODED  BY  R.L.  ARMACOSUMUOOl  10 


0002 

IMPLICIT  RE AL*8 1 A-H.0-2 1 

LMU00160 

0001 

RE AL*6  RHOI N.RATIO.EPS I > THET AO. XEP1 . XE P2 

LMUOOISO 

0006 

COMMON/SHARE /X 1 6SI.0ELI65I,A|6S.65I ,N,M, MN.  NP1 . NM1 

LMU00160 

0005 

COMMON/VALUE/F  *G>  P0.RSIGMAtRJI9OI.RHO 

LMU00170 

0006 

COMMON/E OAL/  H,  Hit  M2 

LMUOOIBO 

0007 

COMMUN/CRST/DEL  XI65I tDELXOI 651 tRHOlNt RATIO. EPSI .THETAO. 
1NTCTR.NUMIN1 ,X1 1 651 .X21 65 1 , XI 1651 . XR21 651 .XR1 1 651 >PR1 , 
2PR2.Pl.FlvRjl|90JfD0TTf PGRA0I65I tDI AGI65I • 
IPREVl.AOELX.RSIGl.Gl.NPHASEcNSATIS 

LNU00190 

LMU00200 

LMU002I0 

LMU00220 

oooe 

DIMENSION  D EL MU 1 901 v DOI 90I.KTESTI65I 

LNU00230 

0009 

MM2  ■  M  ♦  M2 

L  MU 00260 

0010 

GO  TO  (1.2.3,31,  IND 

L  MU00250 

C  INO 

«  0 

LMU00260 

0011 

DO  100  l>l.NM2 

L  MU00270 

0012 

100 

DELMUI 1 1  -  0. 

L MU 00280 

0011 

RETURN 

LNU00290 

C  1ND 

•  1 

LMU00300 

0016 

1 

00  50  I-1.NM2 

LNUOOIIO 

0015 

50 

DELMUI 1 1  ■  RJI 1 1 

LMU00320 

0016 

RETURN 

LMU00330 

C  INO 

«  2 

LMU00160 

0017 

2 

00  60  1*1 , MM2 

LMU00350 

001  B 

60 

DELMUI I i  *  (DELMUI 1 1  -  RJI 1)1 /DEM 

LMU00360 

0019 

RETURN 

LMU00370 

0020 

1 

DO  70  1-1. MM2 

LMUOOIBO 

0021 

CALL  GRADltll 

LMU00390 

0022 

CALL  RESTNT 1 1 , VAL 1 

LNU00600 

0021 

SUM  •  0. 

LMU00610 

0026 

DO  71  JJ-l.N 

LMU00620 

0025 

71 

SUM  -  SUM  ♦  DEL  1 J JIPDELXI JJI 

LMUQ063G 

0026 

1FI INO. ED. 6)  GO  TO  80 

LMU00660 

0027 

OU (II  —(SUM  *  DELMUII 1  l»RH0/IVAl*»2l 

LMU00650 

0028 

GO  TO  70 

LMU00660 

0029 

BO 

DUIII  «  2 *•  (  SUM  ♦  DELMUI  III/T.HO 

LNU00670 

0010 

70 

CONTINUE 

L MOO 0680 

0011 

RETURN 

LNU00690 

0032 

END 

LMU00500 

Figure  21. — Subroutine  LMULT. 
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MAIN 


DATE  •  >02*2 


12/12/30 


C  NAI00060 

C  AUGUST  1971  NAI00070 

C  NAI 00080 

C  MAIN  IS  THE  PROGRAM  THAT  INITIATES  THE  SUNT  ALGORITHM.  THE  INPUT  OF  MAI00090 
C  PARAMETERS.  OPTIONS.  AND  STARTING  POINT  IS  DONE  IN  MAIN.  AFTER  THE  HAIOOIOO 
C  SOLUTION  OF  ONE  NLP  PROBLEM  MAIN  LOCKS  FOR  DATA  FOR  ANOTHER  NLP  PROB.  MAI 001 10 

c  main  has  been  modified  by  armacost'ivtbi  to  include  sensitivity  NA100120 

C  ROUTINES  AMO  BY  GHAEN1I1979I  TO  tNCLUOE  BQUNO  CALCULATIONS  MAI00130 


0001 

IMPLICIT  REAL*8IA-H,0-2I 

MAI  00140 

0002 

REALP4  RHOI N, RATIO, EPSI, THETAO.XE PI, X£P2,TMMAX 

MAI  00150 

0003 

COMMON/ SHARE /X 1 45 1, DEL  1451. A 145. 451 ,N.M ,NN, NPl ,NMl 

MAI00160 

0004 

COMMON  /EOAL/  H,  HI,  M2 

MAI  001 70 

0005 

COMMON  /OPTNS/  NT  1 ,NT2,NT1,NT4,N. 5,NT6,NT7,NT8 , NT9, NT10 

MAI00180 

0006 

COMMON /VALUE/F,G,PO,RS I GMA,RJ(9GI,RH0 

HAI00190 

000’ 

CQMMQN/CRST/DELXI45I .DELX0145 I, RHOI Nt RAT 10, EPSI, THE T AO, 

MAI  00200 

1RSIG1,G1,X1(45I,X2I45I,X1I4SI,XP2I45I.XR1I45I,PR1, 

NAI00210 

2PR2,P1,F1,RJ1190I,OOTT,PGRAO(45I,OI AGI45I , 

NA 100220 

3  PREV3.A0ELX,  NTCTR,  NUMIN1,  NPH4SE,  NSATIS 

MAI  00230 

0008 

COMMON/SEN/ PARI  451, OPAR 145 1.NPAR.I SENS 

MAI  00240 

0009 

COMMON /EX POP T/NE  XOP 1 .NEX0P2, NEX0P3, NEX0P4 ,NEX0PS,XEP1 ,XEP2 

MAI00250 

0010 

C0MMUN/ABG/LV,L2, PER  1451 

MAI  00260 

0011 

COMMON/ ABGl/ FE 1 ,F E2 

MAI  00270 

0012 

COMMON/ ABG  2/ OF  1 1  451  •  0F2 1 45 1 

MAI00280 

0011 

C0NMUN/QP6/NE  X0P6 

MAI 00290 

0014 

COMMON/OP //NE XOP 7 

MAI  00300 

0015 

COMMON/ ABG4/ XV 1 45 1 

MAI00310 

0016 

COMMON/ ABG6/ XX 1 45 1 

MAIO0320 

0017 

DIMENSION  X2I451 

MAI  003 30 

c 

PARAMETER  CARO 

MAI  00340 

0018 

10  REA0(5,50,END«40I  EPSI . RHQIN, THE! AO, RATIO, THMAX.N.N.NZ 

MAI  00350 

0019 

LY  *0. 

MAI00360 

0020 

LZ-O. 

MAI  00370 

0021 

RRft*RHOIN 

NAI00180 

0022 

20  CONTINUE 

MAI  00 390 

0023 

IFILY  .EO.  0.1  GO  TO  IS 

MAI  00400 

0324 

IFINEXQP3  .EQ.  0.1  GO  TO  40 

MAI  00410 

0025 

15  CALL  SET  (TNNAXI 

MAI  00420 

c 

INITIAL  X  VEKCTOR  CARO  FORMAT 

HAI00410 

0026 

IFILY  .GT.  0.1  GO  TO  200 

NA 100440 

0027 

READ  15.601  «X( II .1-1, N) 

MAI  00450 

0028 

00  4000  l-l,N 

MAI  00460 

0029 

XZIl l-XIII 

MAI  00470 

0030 

4000  continue 

MA100480 

0031 

200  RHOIN-RRR 

MAI  00490 

0032 

NTCTR-0 

MAI  00500 

0033 

NP1«N*1 

MAI  00510 

0034 

NM1-N-1 

MAI00520 

c 

SUBROUTINE  REAOIN  IS  UNDER  PROGRAMMER  CONTROL 
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The  vector  S  ,  along  which  the  search  for  the  minimum  is  made,  is  de¬ 
termined  in  subroutine  XMOVE.  The  method  used  to  locate  the  minimum 
along  the  vector  is  the  Golden  Section  search  method,  which  is  based  on 
the  Fibonacci  search  method  (see  [10,  p.  193]). 

The  Golden  Section  method  is  a  one-dimensional  search  method  to 
find  the  minimum  of  the  function  that  does  not  require  the  computation 
of  derivatives.  Figure  23  shows  0PT. 

6.16  0UTPUT 

Subroutine  0UTPUT  contains  most  of  the  "write"  statements  used  to 
print  out  information  on  the  results  of  solving  an  NLP  problem.  It  is 
used  to  print  out  information  after  each  iteration  and  also  to  print  out 
the  solution  estimates  and  the  estimates  of  the  Lagrange  multipliers. 
0UTPUT  is  given  in  Figure  24. 

6.17  PARDIF 

The  subroutine  PARDIF  is  called  by  subroutine  SENS  to  assign  a 
differencing  increment  for  use  in  the  central  differencing  formulas 
indicated  in  Steps  2  and  4  of  Algorithm  2.1.  It  assigns  a  fixed  value 
to  the  differencing  interval  when  the  sensitivity  analysis  is  conducted 
at  the  final  subproblem,  or  assigns  a  value  over  a  specified  interval 
for  sensitivity  analysis  at  the  final  subproblem.  When  a  trajectory 
sensitivity  analysis  is  performed,  PARDIF  returns  a  differencing  inter¬ 
val  which  is  dependent  on  the  particular  subproblem  involved,  refining 
the  value  as  the  subproblem  approaches  the  final  one.  Subroutine  PARDIF 
is  shown  in  Figure  25. 

6 . 18  PERT 

Subroutine  PERT,  called  by  the  program  MAIN,  performs  the  book¬ 
keeping  regarding  various  parametric  changes  and  adjustments  required  in 
the  process  of  bound  calculation.  Figure  26  is  subroutine  PERT. 
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OUTFUT  FRINTS  OUT  INFORMATION  ON  THE  RESULTS  OF  EACH  ITERATION  AND  THEOUTOOIOO 

OUTOOI 10 
OUT 00120 


C 

C  SOLUTION  ESTIMATES  A  NO  THE  ESTIMATES  OF  THE  LAGRANGE  MULTIPLIERS 
C  THIS  ROUTINE  MAS  MOOIFIEO  BY  GHAENII1979I  MHILE  INCORPORATING 
C  BOUND  CALCULATION  TECM1I0UE. 

IMPLICIT  REAL*B(A-H,0-Z> 

REAL *4  RHOIN, RATIO, EPSI.IHETAO 

COMMON /SHARE /X 1 45 1 . DELIAS )»«I45,45) tNtN.MN.NPl.NMl 
COMMON  /EQAL/  H.  HI i  NZ 

CONNUN  /OPTNS/  NT1.NT2.NT3.NT4.NT5.NT6.NTT.NT8.NT9.NT10 
C0MMON/VALUE/F.G.P0.RSIGNA.RJI90I.RH0 

COMMON/CRST/OEL X( 45 I.DELXOI 451. RHOIN. RATIO, EPS  I .THE TA3. 
1RSIG1.G1.RK4SI  .X2I451.  X3I45I  .XR2I45)  .XR1I45I  .PR1 . 

2PR2.Pl .F1,RJ1I90I,OOTT,PGRADI45I.OIAGI45I, 

3  PREV3.A0ELX.  NTCTR.  NUMIN1.  NPHASE.  NS AT  IS 
COMMON/ABG/LV,LZiPER(45 I 
COMMON/ ABGl/FEl .FEZ 
COMMON/ ABG4/ XYI 451 
COMMON/ ABG2/0F1 1451 .0F2I45I 
COMMON/ 0P6/NEX0P6 
C0MMQN/ABG5/TXI 451 .T0ELXI45I 
COMMON/ ABG6/XXI 451 
NZ«M4MZ 

GO  TO  110,201,  X 

10  MR1TE  16,601 

WRITE  16,701  NTCTR. DOTT .RHO.ADELX, NPHASE 
20  WRITE  16,801  F.P0.G.RS1 GM A.H 

IFILY  .GT.  0.)  GO  TO  200 
00  210  lal,N 
XVII l-XII) 

210  CONTINUE 
FE1-F 
GO  TO  220 
200  FE2«F 

220  WRITE  16,901  IXI 11,1-l.N) 

C  CALCULATION  OF  GP  VARIABLES 

IFINEX0P6.NE.1I  GO  TO  315 
INOEX-O. 

CALL  TRANSI INDEX) 

WRITE! 6*  300) (XI l),lal,N) 

00  410  l>l,N 


OUTOOI 30 
OUT00140 
OUTOOI 50 
OUTOOI60 
OUTOOI TO 
OUTOOI SO 
OUT  001 90 
OUT 00200 
OUT 002 10 
0UTO022O 
OUT 002 30 
OUT 00240 
0UT00250 
OUT 00 260 
OUT 002 70 
0UT00280 
OUT00290 
OUT 00 300 
OUTOOI 10 
OUT 00320 
OUT 00330 
OUT 00 340 
OUT 00350 
OUT 00 360 
OUT 00370 
OUT 00  380 
OUT 00390 
QUT00400 
OUT 00410 
OUT  00420 
0UT00430 
0UT00440 
0UT00450 
OUT 00460 
OUT00470 
OUT 004 80 
OUT00490 
0UT00500 
OUTOOSIO 
OUT00520 
OUT 00530 
OUT 00540 
OUT 00550 
OUT 00560 
OUT 00570 
OUT 00580 
OUTOO590 

FORMAT  (50H0*********************P6***P****P*******P***  )  DUT00600 
FORMAT  li0X,6HP0INT«l4.6X,6H  OOTI-E 15.7.6X, 4HRH0-E 15. T,6X, 1OHMAGNIOUX0O61 0 
1TUDE-E15.7.6X.6HPHASEM2I  DUT00620 
FORMAT  <ax,2HF>E15.r,SX,2HP>E15.T,SX,2HG*EI5.7,SX,7HRSI6NA>ElS.7,S0UTOO630 


NT  2 


410  X(1)>TX(I| 

315  WRITE  (6.1101 

GO  TO  130,40). 

30  WRITE  16,1201 

WRITE  16,1001  IRJII),I*1,NZ) 

GO  TO  50 

40  WRITE  (6,100)  I RJI I ) .1*1. NZ) 

300  FORMAT!/, 6X, 'CORRESPONDING  VALUE  OF  GP  VARIABLE  I  S' ./. (6E20.T I ) 
50  RETURN 
C 

60 

70 

80 

90 
100 
110 
120 


1X,2HH>E15.7) 

FORMAT  (6X.25HTHE  CURRENT  VALUE  OF  X  IS/I6E20.7II 
FORMAT  I6E20.T) 

FORMAT  (6X.21HTHE  CONSTRAINT  VALUESI 

FORMAT  I2BX.34HN0T  INCLUOING  THE  NON-NEGATIVITIES) 

ENO 


OUT 00640 
OUT006S0 
OUT00660 
OUT 00670 
OUT 00680 
OUT 00690 


Figure  24. —Subroutine  0UTPUT. 


j 


-  89  - 


T-434 


FORTRAN  IV  C  LEVEL  21 


PARDIF  OATE  ■  80242  12/09/)$ 


0001 


0002 
0003 
0009 
000$ 
0006 
ooo  r 
0008 
0009 
0010 
0011 
0012 


SUBROUTINE  PARDIF 

1$  NARCH  1972 

THIS  SUBRUUT INE  RETURNS  A  DIFFERENCING  INTERVAL  FOR  DPARI20I  WHEN 
CALLED  by  SUBROUTINE  SENS.  ISENS  CONTROLS  THE  VALUE  OF  THE  INTERVAL 
NITH  A  HINIHJM  VALUE  IN  THE  NEIGHBORHOOD  OF  1.E-1S  OE PENDING  ON  KEPI 
ISENS  DEPENDS  ON  THE  KINO  OF  SENSITIVITY  ANALYSIS  BEING  CONDUCTED. 
THE  SUBROUTINE  MAS  CDOEO  BY  R.  L.  ARNACOST. 

IMPLICIT  REAL*8IA-H.0-ZI 
REAL**  KEPI 

COMMON/EXPQPT/NE X0P1 .NEX0P2. NEX0P3. NE XOP4-NEXOP9. XEP1 . XEP2 

CONNQN/SEN/PARI LSI. DP ARI49I.NPAR, ISENS 
MULT  •  NINOI ISENS. 81 
XXEP1  -  XEPl  /  1 1 0.**MULT 1 
IFIXXEPl.LE.O.)  XXEP1  •  l.E-10 
00  10  1-l.NPAR 
10  OPARI 1 laXXEPl 
RETURN 
END 


PAR00040 
PAROOOSO 
PAROOObO 
PAR00070 
PAROOOBO 
PAR 000 90 
.PAR00100 
PAR 001 10 
PAR00120 
PAR00130 
PAR00140 
PAR  00 1  $0 
PAR 00160 
PAR 001  TO 
PAR 00180 
PAR 00190 
PAR 00200 
PAR00210 
PAR 00220 
PAR 00230 
PAR 00240 


Figure  25. —Subroutine  PARDIF. 


FORTRAN  IV  G  LEVEL  21  PGRf  DATE  »  80242  12/13/33 


C  SUBROUTINE  PERT  IS  COOED  BY  GHAENI I 19T9I  TO  PERFORM  THE  PARAMETER  PER00060 

C  AOJUASTNENTS  reouireo  at  various  STAGES  OF  BOUND  CALCULATIONS.  peroooto 


0001 

SUBROUTINE  PERT 

PER00080 

0002 

implicit  aeal»bia-m,o-zi 

PER00090 

0003 

COMHON/ABG/LY. lz. PERILS  1 

PER00100 

0004 

COMMON/SEN/PARI  451 1 DP ARI4SI.NPAR. ISENS 

PER001 10 

000$ 

20 

lz-lzm 

PER00120 

0006 

IF (LZ  .GT.  NPARI  GO  TO  30 

PER00130 

OOOT 

1FIPERILZI  .NE.  0.1  GO  TO  10 

PER00140 

0008 

GO  1(1  20 

PER001S0 

0009 

10 

PAR(LZI-PARILZI»PERILZI 

PER00160 

0010 

30 

RETURN 

PER001 TO 

Qou 

END 

PERO0I80 

Figure  26. — Subroutine  PERT. 


6.19  PEVALU 


Subroutine  PEVALU  is  used  to  compute  the  value  of  the  w  function 
and  the  G  function.  It  makes  use  of  the  values  of  the  objective  func¬ 
tion  and  the  constraint  function,  which  must  have  been  computed  before 
PEVALU  is  called.  The  G  function  is  interpreted  as  the  value  of  the 
dual  objective  function  when  the  W  function  has  been  minimized  for  a 
given  value  of  r  (RH0) .  PEVALU  is  given  as  Figure  27. 
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FORTRAN  IV  6  LEVEL  21  PEVALU  DATE  •  80242  12/11/1* 


0001 

SUBROUTINE  PEVALU 

PEV00040 

c 

PEV00050 

c 

OCTOBER  1970 

PEV00060 

c 

PEVOOOTO 

C  PEVAIU  COMPUTES  THE  VALUE  OF  THE  PENALTY  f UNCTION  AND  THE  VALUE  OF 

THEPEVOOOBO 

C  OUAl  USING  PREVIOUSLY  COMPUTED  VALUES  FOR  F,  AND  RJ. 

PEV00090 

0002 

IMPLICIT  REAL»8IA-H,0-ZI 

PEV00100 

oooi 

REAL*4  RH01 N. RAT IO,EPSI,THETAO 

PEV00110 

0004 

CUMMON/ SHARE  /XI4*I»DEL14SI,AI45,4*1 ,N,M, MN, NP 1 »NN1 

PEVOOI20 

ooos 

COMMON  /EQAL/  H>  HI,  M2 

PEV00110 

0006 

COMMON  /OPTNS/  NT1,NT2,NT},NT4,NTS,NT*,NTTIN!8,NT9,NT10 

PEV00140 

0007 

COMMUN/VALUE /F , G , PO , R S 1 GNA , R J 1 90 1 , R HO 

PEV001S0 

0006 

COMMON/CRST/DELX(4*I,OELXOI4*I,RHOIN,RATIO,EPSI,THETAO, 

PEV00160 

1RSIG1,G1,X1(45I,X2(4S),X1(4SI,XR2(4SJ .XR1I4SI ,PRl, 

PEV00I70 

2PR2,P1,F 1«RJ1I90I,00TT, PGR ADI  451,01 AGI4SI, 

PEV00180 

1  PREV1, ADELX,  NTCTR,  NUMIN1,  NPMASE,  NSATIS 

PEV00190 

0009 

H-0.0 

PEV00200 

0010 

RSI GMA-0.0 

PEV03210 

C  NONNEGS  IF  INCLUDED  ARE  ADDED  TO  P — ARE  POS.  IN  ALL  PHASES 

PFV00220 

0011 

GO  TO  110,101,  NT2 

PEV00210 

0012 

10 

DO  20  1-1, N 

PEV00240 

0011 

20 

RS1GNA«RSIGMA-RH0*DL0G<  XI  III 

PEV00250 

0014 

1G 

GO  TO  140, *0,1*01,  NPHASE 

PEV00260 

C  OBJECTIVE  FUNCTION  -  SIGMA  VIDL.  CL4STS. 

PEV00270 

001* 

40 

F«0.0 

PEV00280 

0016 

SO 

IF  (Ml  100,100,60 

PEV00290 

0017 

69 

00  90  J-l.N 

PEV00100 

0016 

IF  (RJ(Jl)  80,80,70 

PEV001I0 

0019 

70 

RSIGNA«RS1GNA-RHU*DL0G( RJI Jl 1  „ 

PEV0012U 

0020 

GO  TO  90 

PEV00110 

0021 

•0 

F-F-RJIJI 

PEV00140 

0022 

90 

CONTINUE 

PEV00190 

C  EQUALITIES  NOT  AOOED  IN  FEAS.  PHASE 

PEV00140 

0021 

100 

CONT 1  HUE 

PEV00170 

0024 

IF  (M2I  140,140.110 

PEV00180 

002* 

110 

GO  TO  (140,120,1*01,  NPHASE 

PE VOO 190 

0026 

120 

00  110  1-1. MZ 

PEV00400 

0027 

R-N»l 

PEV00410 

0028 

110 

H-H*KJIKI*«2 

PEV00420 

0029 

H-H/RHO 

PEV00410 

0010 

140 

HS-H4R SIGMA 

PE V 00440 

0011 

PO-FPHS 

PEV00490 

0012 

HM S- 2 . PH-RH06DF L 0 AT ( M 1 

PEV00460 

0011 

G-F.HMS 

PEV00470 

0014 

lFINT2.EQ.il  G-G-RHOPOFLOATINI 

PEV00480 

001* 

1*0 

RETURN 

PE  V00490 

0016 

END 

PEV00500 

Figure  27. — Subroutine  PEVALU. 


6.20  PRESEN 


This  subroutine  calculates  the  gradient  of  the  optimal  value 
function,  using  the  gradient  of  the  Lagrangian  when  taken  with  respect 
to  the  parameters  (Step  9  of  Algorithm  2.1).  In  addition,  PRESEN  (when 
invoked)  performs  a  preliminary  screening  of  the  problem  parameters  to 
which  the  optimal  value  function  is  practically  insensitive.  The  cri¬ 
terion  used  for  this  purpose  is  that  if  the  change  in  optimal  value 
function  is  less  than  0.1  percent  of  its  current  value  as  a  result  of 
the  introduced  increment  of  a  given  parameter,  then  that  parameter  is 
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eliminated  from  further  consideration.  This  parameter  screening  is 
invoked  when  the  fifth  variable  in  the  second  option  card  takes  a  value 
of  zero,  as  shown  in  Table  4.  Figure  28  lists  PRESEN. 

6.21  PUNCH 

When  a  call  to  PUNCH  is  made,  the  current  value  of  the  array  x 
is  punched,  along  with  some  of  the  control  cards  that  can  be  used  to 
restart  SENSUMT.  Subroutine  PUNCH  is  shown  in  Figure  29. 

6.22  REJECT 

Subroutine  REJECT  puts  values  of  the  point  at  x  and  the  asso¬ 
ciated  values  of  the  objective  function,  constraint,  functions,  and  W 
function  back  into  their  normal  locations.  These  values  were  temporari¬ 
ly  stored  in  other  locations  by  subroutine  ST0RE.  ST0RE  appears  as 
Figure  30. 

6.23  RH0C0M 

This  subroutine  is  used  to  compute  the  initial  value  of  r  (RH0) 
for  each  NLP  problem.  As  previously  stated,  the  search  for  a  feasible 
starting  point  can  result  in  a  sequence  of  NLP  problems.  The  initial 
value  of  r  is  computed  by  the  formula  specified  for  the  use  of  option 
1  (NT1).  RH0C0M  is  given  as  Figure  31. 

6.24  SEC0ND 

SEC0ND  is  used  to  query  about  the  computer's  clock.  This  is  made 
possible  by  an  internal  clock  which  is  called  by  this  subroutine.  The 
elapsed  time  obtained  and  monitored  by  this  subroutine  is  used  by  the 
subroutines  TIMEC,  TCHECK,  and  SET.  Figure  32  is  subroutine  SEC0ND. 
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FORTRAN  IV  S  LEVEL  21  PRESEN  OATE  •  802*2  12/18/SO 


0001 

C 

SUBROUTINE  PRESENIOU.KTESTI 

Wf  00060 
PRE00070 

c 

c 

1  N6RCH  1976 

PREOOOBO 
PRE 00090 

c 

SUBROUTINE  PRESEN  IS  USED  TO  CALCULATE  THE  GRADIENT  07  THE  OPTIMAL 

PREOOIOO 

c 

value  function  using  the  gradient  of  the  LAGRANGIAN  TAKEN  WITH 

PREOOI 1 0 

c 

RESPECT  TO  THE  PARAMETERS.  PRESEN  IS  CALLED  BV  SENS  MHEN  VARIABLE 

PRE001Z0 

c 

NEX0P5  *0  UR  >1.  ADDITIONALLY.  MHEN  NEXOPS  •  0.  THE  PARAMETERS  HHICHPREOODO 

c 

AFFECT  THE  OPTIMAL  VALUE  FUNCTION  8V  LESS  THAN  0.001  TIMES  ITS 

PRE00140 

c 

CURRENT  VALUE  ARE  ELIMINATED  FROM  FURTHER  CONSIDERATION.  THIS 

PRE00150 

c 

SUBROUTINE  HAS  COOED  BV  R.  L.  ARHACOST. 

PRE0UI60 

0002 

IMPLICIT  REAL*B(A-H,D-2I 

PRE00170 

000) 

COMMON/SHARE/XI 451 .DELIAS)* Al 45 .451 .NtM.MN.NPl.NMl 

PREOOI BO 

0004 

COMMON/VALUE/F .G.P0.RSIGMA.RJ190I.RH0 

PRE 00190 

ooos 

COMMON/EOAL/  H,  HI. M2 

PRE00200 

0008 

COMMON/SEN/ PARI  45 1.DPAR 1451 ,NPAR,I SENS 

PRE00210 

0007 

COMMON/ ABC/ LY . L  2 • PER (45 1 

PRE 00220 

0008 

COMMON/ ABC1/FE1.FE2 

PRE00230 

0009 

C0MM0N/ABG2/DF 11451. DF2I45I 

PRE00240 

0010 

DIMENSION  GXI90I .DU190I »K TEST 1 451 .KL1STI45I 

PRE  00250 

0011 

MPM2  *  N  ♦  M2 

PRE 00260 

0012 

FTEST  •  0.001  •  DABSIFI 

PRE00273 

001) 

DO  100  J-l.NPAR 

PRE00280 

0014 

KTESTIJl  •  0 

PRE00290 

0018 

PARIJI  -  PARI J 1  ♦  OPARIJI 

PREOOIOO 

0016 

CALL  RESTNT I O.OF 1 

PREOO) 10 

0017 

IFIHPM2.EQ.0i  GO  TO  20 

PRE 00)20 

0018 

00  10  1*1 .MPM2 

PRE00))0 

0019 

CALL  RESTNT ( I ,OUI III 

PRE00340 

0020 

10 

CONTINUE 

PRE 00350 

0021 

20 

DEN  -  2.  *  OPARIJI 

PRE  00360 

0022 

PARIJI  •  PARIJI  -  DEN 

PRE00370 

002) 

CALL  RESTNT (0. XFI 

PRE00380 

0024 

IFIMPM2.EQ.0I  GO  TO  40 

PRE00390 

002) 

DO  )0  1*1 *HPN2 

PRE 00400 

0026 

CALL  RESTNT 1 I. GXI III 

PRE00410 

0027 

)0 

CONTINUE 

PRE00420 

0028 

40 

OFEPS  -  IDF  -  XF) /DEN 

PRE00430 

0029 

IFIMPM2.EQ.0I  GO  TO  60 

PRE00440 

00)0 

DO  50  l«l.NPM2 

PRE00450 

00)1 

so 

Dull)  *  (Out  It  -  GXI 1 1 1 /OEM 

PRE00460 

no)2 

60 

SUM  *  OFEPS 

PRE 004 70 

00)) 

1FIH.EQ.0I  GO  TO  80 

PRE 00480 

00)4 

DO  70  1*1 tM 

PRE 00490 

00)5 

70 

SUM  *  SUM  -  RHO/RJIII*DU(II 

PRE 00500 

00)6 

80 

1FIM2.EQ.0I  GO  TO  95 

PRE00S10 

00)7 

TSUN  •  0. 

PRE00520 

00)8 

DO  90  1*1, M2 

PRE 00530 

00)9 

IN  «  1*M 

PRE  00540 

0040 

90 

TSUN  •  TSUN  *  RJI IM)*DU( I  Ml 

PRE00550 

0041 

SUN  •  SUM  ♦  TSUM  •  2./RHQ 

PRE  00560 

0042 

95 

DELI J)  *  SUM 

PRE 00570 

004) 

PARIJI  *  PARIJI  ♦  OPARIJI 

PRE00580 

0044 

OTEST  ■  DABSIDELIJII 

PRE00590 

0045 

1FI DTE  ST. GE. FTEST)  KTESTIJl  *  I 

PRE 00600 

0046 

100 

CONTINUE 

PRE 00610 

0047 

HRITEI6,600I 

PRE00620 

0048 

600 

FORMATI22X. )4H0PT1NAL  VALUE  FUNCTION  SENSITIVITY  //I 

PRE00630 

0049 

00  200  l*l,NPAR,5 

PRE 00640 

0050 

ll*NIN0ll*4,NPAR> 

PRE 00650 

0051 

HR ITE (6,601 1  1 1 JJtDELI JJ) I,  JJ-I.III 

PREO066O 

0052 

601 

FORMAT!  5ITH  DF/DAI  ,  12. 2HI-, G14.7) I 

PRE 00670 

005) 

200 

CONTINUE 

PRE00680 

0054 

IFILY  .GT.  0.)  GO  TO  500 

PRE 00690 

0055 

DO  400  J*1,NPAR 

PRE 00700 

0056 

OF  11 JI*OEL( J) 

PRE 00710 

0057 

400  CONTINUE 

PRE00720 

0058 

GO  TO  700 

PRE00730 

0059 

500  00  610  J*1,NPAR 

PRE 00740 

0060 

Of 21 J)*DELIJI 

PRE 00750 

0061 

610  CONTINUE 

PRE 00760 

0062 

IFILY  .EQ.  0.1  GO  TO  700 

PRE00770 

006) 

PARIL2I*PARIL2I-PERIL2I 

PRE00780 

0064 

RETURN 

PRE  00790 

0065 

700  JJ  *  0 

PREOOBOO 

0066 

00  250  J*1,NPAR 

PRE00810 

0067 

1FIKTESTIJI.EQ.0I  GO  TO  250 

PRE 00820 

0068 

JJ  •  JJ  ♦  1 

PRE00B30 

0069 

KLISTIJJI  ■  J 

PREOO 840 

0070 

250 

CONTINUE 

PRE 00850 

0071 

IFIJJ.EQ.OI  GO  TO  300 

PRE 00860 

0072 

HRITEI 6,602 1 

PRE  008  70 

007; 

602 

FORMAT!  /SIN  DETAILEO  SENSITIVITY  RESULTS  FOLLOH  FOR  PARAMETERS 

IPRE00880 

0074 

HRITEI6,60)I  I K  L I  ST  111,  1*1, JJI 

PRE  00890 

0075 

60) 

FORMAT  I IH  40II2.2H  ,11 

PRE00900 

0076 

HRITEI6.604I 

PRE00910 

0077 

604 

FORMAT!/) 

PRE  00920 

0078 

RETURN 

PRE 009)0 

0079 

>00 

HRITEI6.605I 

PRE00940 

0080 

605 

FORMAT!  42H  THERE  ARE  NO  DETAILED  SENSITIVITY  RESULTS  //) 

PRE 00950 

0081 

RETURN 

PRE  00960 

0082 

END 

PRE00970 

Figure  28. — Subroutine  PRESEN. 
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FUAT KAN  IV  G  LEVEL  21  PUNCH  DATE  •  802*2  12/12/21 


0001 

SUBROUTINE  PUNCH 

PUN00040 

c 

PUN00J50 

c 

OCTOBER  1970 

PUN0QJ60 

c 

PUN00070 

c 

THIS  SUBROUTINE  PUNCHES  THE  STOPPING  POINTS  AND  ASSOCIATED  PARAMETERS 

PUN00083 

c 

SO  THAT  ANOTHER  RUN  MAT  BE  MADE  STARTING  WHERE  THE  CURRENT  ONE 

PUN00090 

c 

STOPPED 

PUN00100 

c 

this  routine  is  called  if  nt**2. 

PUNJOIIO 

0002 

IMPLICIT  REAL*8IA-H,0-ZI 

PUN001 20 

0003 

REAL**  RHOIN.RATIO.EPSI , T HET AO , XEPl ,XEP2 , T 

PUH30130 

000* 

COMMON/ SHARE/X 1 *51 i DEL  1*5  IgA  1*5 .*51 .N.N.MN,  NP 1 ,NM1 

PUN00I40 

0005 

COMMON  /EOAL/  H,  HI.  M2 

PUN00140 

0000 

COMMON  /OPT NS/  NTl.NI2.NT3.NTA.NrS. NT*. NTT. NTS, NT9.NT10 

PJN00160 

000T 

C0MM0N/VALUE/F.G.P0.RSIGMA.RJI90I.RH0 

PUNOOl 70 

0008 

C0MM0N/CRST/DELXI*5I,DEIX0I*5I,RH0IN.RATI J.EPSI .THETAO, 

pjNOoiao 

1RS1G1,G1.X1(*51,X2(*SI ,X3I*SI .XR2IASI .XR1IA5I ,PR1, 

PUNOOl 90 

2PR2.P1,F1,RJ1I90),D0TT,PGRADI*S>,DIAG(*5>. 

PUN00200 

3  PREV3.A0ELX,  NTCTR,  NUNINI.  NPHASE.  NSATIS 

PUN002 1 0 

0009 

COMMUN/EXPOPI/NEXOPl .NEX0P2.NE X0P3, NE XOP* .NEX0P5 , XEP1 . XE P2 

PUN00220 

0010 

7*60.0 

PUN 002 30 

0011 

WRITE  17.101  EPSI.RHU. THETAO. RATIO. T.M.N.MJ 

PUN00240 

c 

TMMAX.1S  SET  TU  60.  SECONDS 

PUN00240 

0012 

WRITE  17.201  (XIII, 1*1, N> 

PUN00260 

0013 

NT1*3 

PUN00270 

c 

SET  MHO  OPTION  SO  THIS  VALUE  OF  RHO  MILL  BE  USE  FOR  THE  RESTART. 

PUN00280 

001* 

WRITE  17.301  NT1.Nt2.NT3.NTA»NT5.NTA.NT7.NT8»NT9.NTIO 

PUN00290 

0015 

WRITE  17,201  XE  PI ,  XEP2 

PUN00300 

0016 

WHITE  (7.301  NEXOP1.NEXOF2.NEXOP3 

PUN00310 

001/ 

RETURN 

PUNOO  120 

c 

PUNOO 3 30 

0018 

10 

FORMAT  I 5E 12. 5. 31*1 

PUN JO  340 

0019 

20 

FORMAT  I6E12.3) 

PUNOO 350 

0020 

30 

FORMAT  110(71 

PUN00160 

0021 

ENO 

PUN00370 

Figure  29. — Subroutine  PUNCH. 
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0001 

c 

SUBROUTINE  REJECT 

RE JOOOAO 
REJ00050 

c 

c 

OCTOBER  1970 

RE JOOOAO 
REJOOOTO 

c 

RtJECT  RET  URNS  THE  STUREO  VALUES  OF  THE  OBJECTIVE  FUNCTION.  THE 

REJ00080 

c 

CONSTRAINT  FUNCTIONS  ANO  THE  PENALTT  FUNCTION  TO  THEIR  NORMAL  LOCAT IONRE J00090 

0002 

IMPLICIT  REAL*8(A-H.0-ZI 

REJOOIOO 

0003 

REAL**  RHOIN.RATIO.EPSI. THETAO 

REJOOIIO 

0004 

COMMON/SHARE /X  1*5 1.DEL  1  A5I.AI  AS. *51. N.N.MN.’IPl.NMl 

RFJ00I20 

0004 

COMMON  /EUAL/  H,  HI,  M2 

REJ00I30 

0006 

CQMMON/VALUE/F, G.PO.RSIGMA.RJ 1901, RHO 

RE J001A0 

0007 

COMMON/CRST/OELXIA5I ,DE LXOI *5 1 . RHOI N, RATIO, EPSI .THETAO, 
1RSIGI.G1.X1IASI ,X2I*5I,X3IA51,XR2IASI .XR1IA5I ,PR1, 

2PR2.P1 ,F1,RJ1(90I,OOTT,PGRADIASI,OIAGIA5I, 

3  PREVi.ADELX.  NTCTR,  NUMINI,  NPHASE,  NSATIS 

REJ00150 

REJ00160 

REJ00170 

REJOOISO 

0008 

OU  10  1*1, N 

REJ00190 

0009 

10 

XI1I*X1III 

RE J00200 

0010 

MM2*M*M2 

REJ002IO 

0011 

00  20  J«1,MM2 

REJ00220 

0012 

20 

RJIJI-RJ1IJ1 

REJ00230 

001 3 

P0*P1 

REJ002AO 

0014 

RSIGMA*RSIG1 

RE J002S0 

0014 

G-Gl 

REJ00260 

0016 

F*F  l 

REJ002T0 

0017 

H»H1 

RE J002A0 

0018 

RETURN 

RE J00290 

0019 

ENO 

REJ00300 

Figure  30. — Subroutine  REJECT. 
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0001 

SUBROUTINE  RHOCOM 

AH000040 

C 

RH000050 

c 

OCTOBER  1970 

AH000069 

c 

AH090070 

c 

SUBROUTINE  TO  COMPUTE  INITIAL  RHO  VALUE 

AH000080 

c 

CONTROLLED  BY  COL.  7  ON  OPTION  CARO 

RHOOO 190 

0002 

IMPLICIT  REAL68IA-H.0-Z) 

AHOOOIOO 

000) 

REAL *4  RH01N, RATIO,  EPSI  .THETAO 

AH000I10 

0004 

COMMON/ SHARE /X< 45), DEL <451, A (45, 451  , N, M , MN, NP l , NM1 

AHQ30120 

0005 

COMMON  /OPTNS/  N71,NT2,NI),N14,NT5,NT6,NT7,NT8,NT9,NTI0 

RHOOOl 10 

0006 

COMMON/ VALUE/F,  G, PO , AS l GMA,R) ( 901 .RHO 

KH000140 

0007 

COMMON/CRS  7/OEL  X( 451 ,OELXOI4S>, RHOIN, RATIO, EPSI , THETAO, 

RHOOO ISO 

IRSIG1,G1,X1(45I,X2<451.X)<45) ,XR2< 451 ,XR 1(45), PR1. 

RHOOOl 60 

2PR2,P1.F1,RJ1I90I ,OOT7,PGRAOI 451,01 AG! 451, 

RHOOOl 70 

3  PREV3,  ADEL  X ,  NTCTR,  NUMIN1,  NPHASE,  NSATIS 

RHOOOl SO 

0008 

GO  TO  (110,50,10,190),  NT  1 

RHOOOl 90 

0009 

10 

AH0*AH01N 

RHQ00200 

0010 

20 

IF  (RHO)  30,30,40 

AH0002 10 

0011 

30 

RHO* 1 • 

RH000220 

0012 

40 

RETURN 

RH000230 

001) 

50 

NPAR 1*1 

RM000240 

0014 

60 

RHO*  1  • 

RMQ00250 

c 

2  MEANS  RHO  WHICH  MINIMIZES  GRADIENT  MAG. 

RHQ00260 

0015 

CALL  GRAD  121 

RMU00279 

0016 

00  70  1*1, N 

RM000280 

0017 

70 

PGR AO ( 1I*0ELX0I1) 

RH000290 

0018 

RHO* 2. 

RH030309 

0019 

CALL  GRAO  (21 

RHOOO) 10 

0020 

DO  80  l-l. N 

RHO 00320 

0021 

DELXOI II-DELXOI II -PGR ADI 11 

RH000330 

0022 

80 

PGRAOl  ll-PGRADI  1 1-DELXOd  1 

AH000I40 

002) 

GO  TO  (90,1301,  NPAR1 

RHQ00350 

0024 

90 

0071*0. 

RHOOO 360 

0025 

0072*0. 

RH000370 

0026 

00  100  l-l. N 

RH000380 

0027 

0071-O0T14DELX0I I >*PGRADt II 

RHQ00393 

0028 

100 

OQ72-OOT2»DELXOI I 1**2 

RM000400 

0029 

RHO*OABS( DOT  1/00721 

RH000410 

0030 

GO  70  20 

RM000420 

C  3 

MEANS  COMPure  RHO  SO  AS  TO  MINIMIZE  DEL  PI/DOP/I.IDEL  P 

RH000430 

00)1 

110 

NPAR2*L 

RM000440 

0032 

120 

NPAR 1*2 

RH000450 

C 

USE  OF  ANO  DR  SUBROUTINE 

RH000460 

0033 

GO  TO  60 

AH000470 

0034 

130 

RHO*  l  • 

RM000480 

C 

ASSUME  SIGMA  TERN  IS  CONSID.  GRTER  THAN  F  TERN 

RM000490 

0035 

CALL  SECORO  121 

RH000600 

0036 

00  140  1*1, N 

RH000510 

0037 

140 

OELXI 1 >*PGAADI 1 1 

RH000520 

00)8 

CALL  INVERS  111 

RH000530 

00)9 

DO  150  1*1, M 

RH000540 

0040 

XKIfOELXIII 

RHO 00550 

0041 

150 

OELXI I 1*081X0(11 

RH000560 

0042 

CALL  SECORO  121 

RH000570 

0043 

CALL  INVERS  111 

RHOOO 5 00 

0044 

DO  160  1*1, N 

RM000590 

0045 

160 

XR2III*0ELX(I) 

RH000600 

0047 

170 

0071*0. 

RH000620 

0048 

DOT 2*0. 

RH000630 

0049 

DO  180  1*1, N 

RH000640 

0050 

D0T1*00T14PGRADII )*X1I 1 1 

AM000650 

0051 

180 

OOT2*DOT24DELXOI 1 l*XR2l  II 

RH000660 

0052 

RHO* DSQRTI DABSI DOT 1/OOT 2 1 1 

RH000670 

0053 

GO  TO  20 

RH000680 

0054 

190 

NPAR2*2 

RH000690 

C 

RHO  MINIMIZES  2ND  ORDER  MOVE 

RH000700 

0055 

GO  TO  120 

RH000710 

t 

USES  INTERNAL  SUB.  TO  COM  /DDP/-1  DF  ANO  /OOP/-  DR 

RH000720 

0056 

200 

0071*0.0 

RH000730 

0057 

0072-0.0 

RHO00740 

0058 

00  210  1*1, N 

RH000750 

0059 

00T1*X1| 1 |P*2*00T1 

RH000760 

0060 

210 

D0T2-X1I I )*XR2I I 1 *0072 

RH000770 

0061 

RHO*DABS(DOTl/OOT2> 

RH000780 

0062 

GO  TO  20 

RHOOO 790 

006) 

ENO 

RHOOO 800 

Figure  31. — Subroutine  RH0C0M 
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0001 

SUBROUTINE  SECONOISECSI 

SEC000S0 

0002 

CALL  CLOCK! 11 

SEC 00060 

0003 

SECS  •  1 

SEC00070 

000* 

SECS  «  SECS  /  100.0 

SEC00080 

00  os 

RETURN 

SEC00Q90 

ooo* 

ENO 

SEC00100 

Figure  32. — Subroutine  SEC0ND. 


6.25  SEC0RD 

Subroutine  SEC0RD  evaluates  the  matrix  or  second  partials  of  the 
W  function.  It  calls  on  the  user-supplied  subroutine  MATRIX  to  evaluate 
the  upper  triangle  of  the  matrix  of  the  second  partials  of  the  objective 
and  the  constraint  functions.  Subroutine  SEC0RU  is  listed  in  Figure  33. 

6.26  SENS 

This  subroutine  calculates  the  solution  point  and  optimal  value 
function  sensitivities  for  each  subproblem  (if  required),  as  well  as  the 
final  problem,  by  implementing  Steps  1-3  and  7-8  of  Algorithm  2.1  of 
Section  2.  The  various  choices  for  sensitivity  calculations  are  given 
in  options  3,  4,  and  5  of  the  second  option  card  (see  Table  4).  SENS  is 
called  by  the  program  MAIN  and  the  subroutine  B0l)Y.  It  calls  subroutine 
LMUI.T  to  calculate  the  Lagrange  multiplier  sensitivities.  SENS  is  shown 
as  Figure  34. 

6.27  SET 

Subroutine  SET  is  called  once  by  MAIN  at  the  start  of  the  attempt 
to  solve  a  problem.  It  obtains  and  stores  the  value  of  the  computer's 
clock.  It  obtains  the  value  of  the  clock  by  calling  a  system  library  or 
user  coded  subroutine  that  queries  the  computer's  clock  and  returns  its 
value  as  a  floating  point  number  in  seconds.  Figure  35  lists  subroutine 
SET. 
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F0RTA4N 

IV  &  LEVEL  21  SECORO  DATE  -  80242  12/46/12  1 

0001 

SUBROUTINE  SECORO  1 I S  B 

SEC00040 

C 

SEC00050 

C 

OCTOBER  1970 

SEC00060 

c 

SEC00070  j 

C  SECOKI)  EVALUATES  THE  MATRIX  OF  SECOND  PARTIAIS  OF  THE  PENALTY 

SEC00080 

C  FUNCTION. 

SEC 00090 

1 

C 

secooioo 

1 

c- 

(11  MEANS  DONT  COMPUTE  GRAD.  OUTER  PROOUCTMN  SECOROI 

SEC  001 1 0 

0002 

IMPLICIT  R£AL*81A-H,0-2I 

SEC00120 

1 

OOOi 

REAL**  RHOIN, RATIO. EPSI.THETAO 

SEC00130 

1 

0004 

COMMUN/SHAAE/X (951. DEL ( 95 ) . A( 95.951 .N.M.MN.NPl ,NM1 

SEC00140 

1 

0005 

CONNON  /EQAL /  H,  HI,  HZ 

SEC001 50 

1 

0006 

COMMON  /OPTNS/  NTI,NT2,NT3,NT9,NTS,NT6,NT7,NTB,NT9,NT10 

SEC00160 

[ 

0007 

COMMON/VALUE/F.G.P0.RSIGMA.RJI9OI.RHO 

SEC00170 

oooa 

COMM JN/CRST/ DEL X( 951 .OELXOI 951, RHOIN.RATIO.EPS! .THETAO, 

SEC001 80 

1RSIG1.G1 , XI (951 ,X2(95I ,X3(95) . XR2 (95 ( , XR1 (95) ,PR1 , 

SEC00I9O 

! 

2PR2  »Pl»FitRJl(90}»O0TT»  PGR  ADI  451,01  AG  I  451 » 

SEC00200 

3  PREV3. ADELX,  NTCTR,  NUMIN1  NPHASE,  NSATIS 

SEC00210  ! 

i 

0009 

DO  10  1*1 , N 

SEC00220 

00 10 

DO  10  J*I.N 

SEC  002  30 

0011 

10 

A( 1. Jl-O. 

SEC  00240 

0012 

GO  TO  (230.201,  IS 

SEC00250 

CGRAD.  TERM  NOT  PREV-  CONFUTED 

SEC00260  1 

0013 

20 

00  30  1*1, N 

SEC00270 

i 

0014 

DO  30  J*l,l 

SEC00280 

! 

0015 

A (  I. Jl*0.0 

SEC00290 

0016 

30 

CONTINUE 

SEC00300 

0017 

GO  TO  (90,601,  NT 2 

SEC003I0  i 

0018 

40 

DO  50  1*1, N 

SEC00320 

0019 

50 

A(  l.il*RHU/X(ll**2 

SEC00330 

0020 

60 

continue 

SEC00340 

1 

0021 

IF  1 N.LE.OI  GO  TO  130 

SEC00350  \ 

0022 

00  120  IN*1,M 

SEC00360 

0023 

IF  IRJIINI)  120,120,70 

SEC00370  , 

0  024 

70 

CALL  GRAD1  (IN) 

SEC  00380 

0025 

TT  *  RHO/RJI 1N)**2 

SEC00390 

0026 

00  110  1*1, N 

SEC00900 

i 

002  7 

IF  ( DEL (III  BO, 110,80 

SEC00910 

oo2» 

80 

T*TT  *DEL ( 1 1 

SEC00920 

i 

0029 

00  100  J*lfl 

SEC00930 

0030 

IF  (OEL(Jl)  90.100,90 

SEC00990 

0031 

90 

All, J)*A<I.JI»T*DEL(JI 

SFC00950 

, 

0032 

100 

CONTINUE 

SEC00960 

0033 

110 

continue 

SEC00970 

1 

. 

00  39 

120 

CONTINUE 

SEC 00980 

' 

C  EQUALITY  CONSTRAINTS 

SEC 00990 

003a 

130 

IF  (Mil  210,210,190 

SEC00500 

0036 

140 

GO  TO  (210,150,230),  NPHASE 

SEC00510 

0037 

150 

RQ-2./RH0 

SEC00520 

00)8 

DO  200  JU-ltN Z 

SEC00530 

00  39 

IN*N»JJ 

SEC00590 

• 

0040 

CALL  GRA01  (IN) 

SEC00550 

0041 

DO  190  I-l.N 

SEC00560 

! 

0042 

IF  (DELHI)  160,190,160 

SEC00570 

• 

0043 

160 

T*RQ*0ELIII 

SEC00580 

00*9 

DO  130  3*1,1 

SEC00590 

0045 

IF  (OEL(J)I  170,180,170 

SEC00600 

0046 

170 

All , J)*A( I, J)*T*DELI Jl 

SEC00610 

00*7 

ISO 

CONTINUE 

SEC 00620 

0048 

190 

CONTINUE 

SEC 006 10 

0049 

200 

CONTINUE 

SEC00640 

1 

0050 

210 

DO  220  1*1, N 

SEC 00650 

0061 

DIAGI II*A( 1,11 

SEC00660 

0052 

220 

A(  I  .  I  )*0. 

SEC00670 

! 

C  READY  NOR  FOR  MATRIX  OF  2ND  PART  I ALS  OF  RESTRAINTS 

SEC00680 

0053 

230 

GO  TO  (290,510,5201,  NT  10 

SEC 00690 

■  . 

0054 

290 

IF  (M.LE.OI  GO  TO  390 

SECOOTOO 

, 

0055 

DO  330  IN*1,M 

SEC00710 

1 

0056 

LORN* 2 

SEC00720 

! 

C  CONSTRAINT  assumeo  nonlinear 

SEC00730 

f 

0057 

CALL  MATRIX  (IN, LORN) 

SEC 00790 

i 

OOf  8 

IF  ILORN.LT. 21  GO  TO  330 

SEC00750 

f 

0059 

IF  IRJIIN)  ,GT.  0.01  GO  TO  280 

SEC00760 

C060 

00  260  1*2, N 

SEC00T70 

0061 

I M l « 1—1 

SECOOTOO 

0062 

00  260  0*1,  IM1 

SEC00790 

006J 

IF  (AIJ.II)  250,260,250 

SEC00800 

0064 

250 

A(l,JI*Atl,J>-A(J,l> 

SEC 008 10 

0065 

AIJ.II *0. 

SEC00820 

' 

0066 

260 

CONTINUE 

SEC  00830 

006  7 

DO  270  1*1. N 

SFC00840  j 

i 

0068 

DIAGII)*DIAG())-A(1,II 

SEC  00850 

*> 

0069 

2T0 

A(  1,11*0.0 

SEC00860 

0070 

GO  TO  330 

SEC00870 

0071 

2S0 

T— RHO/RJI  INI 

SEC  00880  | 

m 

0072 

00  300  1*2. N 

SEC  00890 

iJMllltl 

0073 

INI* 1“ 1 

SEC 00900 

l 

Figure  33. --Subroutine  SEC0RD. 

J 

Wiu^1  i  u •—» * '.■  '»** '  rr~*'*+ - ^ -**«•»-  -  -v-* ^  "•*■'*  ■*  *  '  ■'■*-i*FW*  "  I  ju ^ 


oo  n 

00  300  J* 1 ,  1  Ml 

SEC00910 

0075 

IF  lAU.III  290.100.290 

SEC 00920 

0076 

290 

A(l,JI*AII,J|4T*AtJ,l) 

SEC00910 

0077 

At J.II-O. 

SEC 00940 

0078 

300 

CONTINUE 

SEC009S0 

0079 

00  320  1*1. N 

SEC 00960 

0080 

IF  IAII.111  310.320.310 

SEC 00970 

0081 

310 

01 AG( 1 1*01 AGi ll«T*A 11.11 

SEC004B0 

0082 

All, 11*0. 

SEC00990 

308  J 

320 

CONI INUE 

SECOIOOO 

0084 

330 

CONTINUE 

SEC01010 

0085 

390 

CONTINUE 

SEC01020 

0086 

GO  TO  1520,150.520).  NPHASE 

SEC0I030 

0087 

350 

IF  1M2.E0.0)  GO  TO  420 

SEC01040 

C — 

EQUAL  1  TV  SECOND  PARTI  ALS  HERE 

SEC01050 

0088 

IF  INT10.GE.2I  CO  TO  420 

SEC01060 

0089 

00  410  11*1. M2 

SEC01070 

0090 

IN*M»II 

SECOIOAO 

0091 

LOAN-2 

SEC01090 

0092 

CALL  MATRIX  (IN, LORN) 

SEC01100 

0095 

IF  ILORN.LT. 21  SO  TO  410 

SEC01110 

0094 

T*2.*RJ( INI/RHO 

SEC01120 

0095 

00  380  1-2. N 

SEC0U10 

0096 

IM 1*  I- 1 

SEC01 140 

0097 

00  170  J-1.IM1 

SEC0I150 

0098 

IF  (AIJ.III  360,370.360 

SEC0I160 

0099 

390 

A(I,JI-AU,J)4T*A(J,!I 

SEC01 170 

0100 

At J, 11*0.0 

SEC  Jl  180 

0101 

370 

CONTINUE 

SEC01190 

0102 

360 

CONTINUE 

SEC01200 

0103 

00  400  1*1, N 

SEC01210 

0104 

IF  (All, 1)1  390.400,190 

SEC01220 

0105 

390 

01 AGI I 1*01 AGI 1 »»T*AI 1,1 1 

SEC01230 

0106 

All, 11*0.0 

SEC01240 

0107 

900 

CONTINUE 

SEC01250 

0108 

910 

CONTINUE 

SEC01260 

C  GET  MATRIX  OF  2ND  PART  SALS  OF  OBJECTIVE  FUNCTION 

SEC01270 

0109 

920 

LLL*2 

SEC01280 

0110 

CALL  MATRIX  (O.LLL) 

SEC01290 

0111 

IF  (LLL.LT.2I  GO  TO  490 

SEC01300 

0112 

00  440  1*2, N 

SEC01110 

0111 

IMl-l-1 

SEC01320 

0114 

00  440  J*1 • IM1 

SEC01310 

0115 

IF  (AIJ.III  430.440.430 

SEC01340 

0116 

930 

A(1,JI-A(I,JI»A(J,II 

SEC01 350 

OUT 

990 

At  J,13»Ml,Jt 

SEC01360 

0118 

00  470  1*1, N 

SEC01170 

0119 

IF  1 A I  I. 1 1 1  450.460,450 

SEC01180 

0120 

950 

A( 1 , 1 1*01 AGI 1)4 A( 1,1) 

SEC01390 

0121 

GO  TO  470 

SEC01400 

0122 

960 

A 1 1 , I l-OI AGI 1 1 

SEC01410 

0123 

970 

CONTINUE 

SEC01420 

0124 

960 

RETURN 

SEC01410 

0125 

990 

00  500  1*1, N 

SEC0I440 

0126 

All, 11*01 AGI II 

SEC01450 

0127 

00  500  J*I.N 

SEC01460 

0128 

500 

AII,JI«A|J,I) 

SEC01470 

0129 

GO  TO  480 

SEC01480 

0130 

510 

GO  TO  1520,350,3501,  NPHASE 

SEC01490 

0131 

520 

00  510  1*2, N 

SEC01500 

0132 

1 191-1-1 

SEC0I5I0 

0133 

DU  530  J-'1,IN1 

SEC01520 

0134 

530 

A  I  J  ,  I  I  •  A 1 1 ,  J  1 

SEC01510 

0135 

00  540  1*1. N 

SEC0I540 

0136 

590 

Al 1 , 1 1*01 AGI 1 I 

SEC01550 

0137 

GO  TO  480 

SEC 01 560 

0138 

ENO 

SEC01570 

Figure  33. — continued 


T-434 


FORTRAN  IV  C  LEVEL  21 


-  99  - 

SENS 
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0001 


0002 

0003 

0004 

0005 

0006 

0007 

0008 


0009 

ooio 

oan 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 

0021 

0022 

0023 

0024 

0025 


0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0043 

0046 

0047 

0048 

00*9 


SUBROUTINE  SENS 


SEN00060 

C  SENOOOTO 

C  1  MARCH  1976  SEN00080 

C  SEN00090 

C  THIS  VERSION  OF  THE  SENSITIVITY  ANALYSIS  SUBROUTINE  IS  USED  TO  SEN00100 

C  COMPUTE  THE  DIRECTIONAL  DERIVATIVES  QF  X  AND  F  HUH  RESPECT  TO  SFNOOUO 

C  CERTAIN  PARAMETERS  CODED  IN  THE  ARRAY  PARI20I.  THE  DIRECTIONAL  SEN00I20 

C  OERIVAI IVES  ARE  ESTIMATED  FUR  ONE  PARAMETER  AT  A  TIME  WITH  NPAR  BE  I NGSENOOl SO 
C  THE  NUMBER  OF  PARAMETERS  INVOLVED  I..'  THE  SENSITIVITV  ANALYSIS.  THE  SENO0I40 
C  USE  OF  THE  PARAMETERS  PARI20I  MUST  A£  CONSISTENT  THROUGHOUT  THE  SFN00I50 

C  USER'S  SUBROUTINES.  SFNOOUO 

C  the  subroutine  is  used  for  a  sensitivity  analysis  at  the  final  sub-  senooito 

C  PROBLEM  OR  FOR  A  SENSITIVITV  ANALYSIS  AT  EACH  SUBPROBLEM  ALONG  THE  SFN00180 

C  MINIMIZING  TRAJECTORY.  0PARI20)  IS  THE  ARRAY  OF  OIFFERENCING  SEN00190 

C  INTERVALS  CORRESPONDING  TO  THE  PARAMETERS  PARI20I.  0PARI20T  IS  SEN00200 

C  ASSIGNED  VALUES  IN  SUBROUTINE  PAROIF.  SEN002I0 

C  THIS  APPROACH  TO  SENSITIVITY  ANALYSIS  IS  OUE  TO  A.  V.  FIACCO.  THF  SEN00220 
C  FIRST  VERSION  HAS  CODED  BY  B.  CAUSEY.  THE  SECOND  VERSION  MAS  COOED  SEN00230 

C  BY  M.  C.  MYLANOER.  THE  THIRD  VERSION  HAS  AN  EXTENSION  OF  THE  SEN00240 

C  SECOND  VERSION  TO  PERMIT  SENSITIVITY  ANALYSIS  ALONG  THE  MINIMIZING  SEN00250 
C  TRAJECTORY,  AND  HAS  COOEO  BY  R.  L.  ARMACOST.  THIS  IS  THE  FOURTH  SEN00260 
C  VERSION  UHICH  INCORPORATES  THE  OPTION  TO  CALCULATE  THE  PARTIAL  SEN00270 

C  DERIVATIVES  OF  THE  LAGRANGE  MULTIPLIERS  {SUBROUTINE  LMULTI  ANO  TO  SEN00280 
C  CONDUCT  A  PRELIMINARY  SCREENING  OF  THE  PARAMETERS  BY  ESTIMATING  THE  SEN00290 
C  FIRST  ORDER  SENSITIVITY  OF  THE  OPTIMAL  VALUE  FUNCTION  ISUBROUTINE  SFN00300 
C  PRESEN).  THIS  VERSION  HAS  COOEO  BY  R.  L.  ARMACOST.  SEN00310 

C  THIS  ROUTINE  TOGETHER  HITH  OTHER  SENSITIVITY  ROUTINE  MERE  SEN00320 

C  REDIMENS IONED  TO  45  I.E.  XI45)  ANO  PARI45I.  FURTHERMORE  SEN00330 

C  IT  MAS  SLIGHTLY  MODIFIED  TO  INTERFACE  BOUND  CALCULATION  SEN00340 

C  ROUTINES.  GHAEMII19T9).  SEN00350 

IMPLICIT  REAL*BI A-HtO-Z )  SENJ0360 

RE AL*4  RHOIN, RATIO, EPSI .THETA0.XEP1.XEP2  SEN00370 

COMMON/ SHARE /XI 45 1 ,0ELI45) » Al 45,45) , N«M,MNt NPltNMl  SEN00380 

COMHON/E  QAL/H, HI , HZ  SFN30390 

CONNON/OPTNS/NT1,NT2,NT3,NT4,NT5,NT6,NT7,NTB,NT9,NT10  SEN00400 

COMMON/ VAL  U£/F»G,PO,RSI GMA.RJI90I «RHO  SEN00410 

COMMON/CRST/OELXI45 I ,0ELXO( 45 ), RHOIN, RAT 10, EPSI .THETAO,  SEN00420 

1NTCTR,NUMINI,X1I4S),X2I45I,X3I45),XR2I45),XR1I45I.PR1,  SEN30430 

2PR2,P1,F1,RJ1(90) ,00TT,PGRA0I4S>,DIAGI45I,  SEN00440 

3PREV3, AOELX.RSIGl.GliNP HA SE,N SATIS  SEN00450 

COMMON/ SEN/ PAR  1451* OP AR 1451 , NPAR, I  SENS  SEN00460 

COMMON/ EXP0PT/NEX0P1 ,NEXOP2,KEXOP3iNEXOP4,NEXOP5,XEPl,XEP2  SEN00470 

C0MM0N/ABG/LY,LZ,PERI45 I  SEN00480 

C0NM0N/ABG1/FE1 >FE2  SEN00490 

COMMON/ ABG2/ OF  1 1 451 , 0F2 1 45 1  SE  N00500 

C0MM0N/ABG5/TXI 43) .TOELXI  45)  SEN00510 

C  0MM0N/0P6/NE  X0P6  SEN00520 

DIMENSION  DELT 145 ) , DELT XI 451 1  X2HI 431  •  X3HI  45)  SEN00530 

DIMENSION  0ELNUI90I t  OUi 901 ,KTESTI 451  SEN00540 

00  5  1*1 , N  SEN00550 

OELTXi tl  •  DELXIII  SEN00560 

5  DELT III  «  OELXOI I )  SEN00570 

CALL  STORE  SEN00580 

AVAL  •  I.0E-10  SEN00590 

MPMZ  *  M  ♦  MZ  SEN00600 

I  SENS  •  I  SENS  ♦  l  SEN00610 

CALL  PAROIF  SEN00620 

C  MRITE  OUT  X,  RHO  ANO  PAR.  SFN00630 

HRITEI 6,101  RHO  SEN00640 

10  FORNATIIH  //30X#20HSENS!T I VI TV  ANALYSIS  //  SEN00650 

15X,22HTHE  VALUE  OF  RIRHOI  IS  ,E12.5I  SEN00660 

MRITEI6.20I  SEN00670 

20  F0RMATI/6X,  64HTHE  POINT  AT  MHICH  THE  ESTIMATE  OF  SENSITIVITY  MILLSEN00680 


1BE  MAOE  IS  I 
00  40  1*1, N,6 
II*M1N0II»5,N> 

MRITEI6.30I  (IJ.XIJI),  J«Z,III 
30  FORMAT  I  6I2X,2HXI ,I2,2H)*,G14.T)  ) 

40  CONTINUE 

CALCULATION  OF  GP  VARIABLE  VALUE 
IFINEX0P6.NE.ll  GO  TO  660 
INOEX«0. 

CALL  TRANSI  INDEX) 

MR1TEI4.400I 

400  F0RMA7I16X, 'CORRESPONDING  POINT  FOR  GP  PROBLEM  IS' I 
DO  440  1*1, N, 6 
II*NIN0II»5,NI 

640  MRITEI6ilOMIJ«XIJI»tJ*l,ll> 

DO  650  I* 1 t N 
450  XI I l*TXI 1 1 
640  CONTINUE 

MRITEI6, 501  II ItPARIl I , 0PAR1 1 II,  I-l.NPAR) 

50  FORMAT I45H0  PARAMETER  VALUE  DIFFERENCING  INTERVAL  / 
1  I14.2X, 614.6, 6X.G14. 51  I 
HRITEI4.33) 

55  FORMAT  I //I 
EVALUATE  F,  G  ANO  H. 


SEN0069Q 

SENOOTOO 

SEN00T10 

SEN00720 

SFN00730 

SEN00740 

SEN00750 

SEN00760 

SEN00770 

SEN00780 

SEN00T90 

SEN00800 

SEN00810 

SEM00820 

SEN00830 

SFN00840 

SEN00850 

SEN00860 

SENOOSTO 

SENOOSBO 

SEN00890 

SFM30903 

SEN009IO 

SEN00920 


Figure  34. — Subroutine  SENS 
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0J50 

CALL  KESTNIIO.F 1 

SEN00930 

0051 

IFINPM2.Ed.OI  GO  TO  85 

SEN00940 

0052 

00  80  J»  1  »  MPM2 

SEN00950 

005 3 

CALL  RESTNT 1 J.RJI 31 1 

SEN00960 

0054 

80  CONTINUE 

SEN00970 

0054 

IFINEX0P5.EQ.21  GO  TO  85 

SEN00480 

0056 

CALL  PRE  SENIOUtKTESTI 

SEN00990 

0057 

IFILY.CT.OI  GO  TO  400 

SENOIOOO 

C  COMPUTE  DEL  F. 

SEN01010 

0058 

85  CALL  GJU01I0! 

SEN01023 

0054 

00  90  iHtN 

SEN01330 

0060 

X1HI  II  •  DELHI 

SENOI 040 

0061 

90  CONTINUE 

SEN0105D 

C  COMPUTE  IOEL 1**2  P  -  STORED  IN  A. 

SEN01060 

00 62 

CALL  SECOAOI 21 

SEN0IO70 

C  PERFORM  THE  L-U  DECOMPOSITION  OF  A. 

SEN01090 

0063 

DO  100  l«l.N 

SFN3I090 

0064 

OELXIi 1-0.0 

SENOIIOO 

0065 

100  CONTINUE 

SEN01110 

0066 

NP-NT3 

SFN01 120 

006  7 

NT  3*  1 

SENOI  UO 

0068 

CALL  INVERSill 

SEN01 140 

C  CHECK  TO  MAKE  SURE  AN  ORTHOGONAL  MOVE 

IS  NOT  ATTEMPTED. 

SENOI 150 

0064 

DO  110  l«l.N 

SEN01160 

0070 

IFIDELXI ll.EQ.O.OI  GO  TO  110 

SENOI 170 

0071 

WRITE <6. 1051 

SENOI 180 

0072 

105  FORMAT (94H0  THE  MATRIX  OF  SEC3ND  PARTI ALS  IS  NOT 

POSITIVE 

OEFINI TESFNOl  190 

It  SENSITIVITY  ANALYSIS  IS  TERMINATED  1 

SEN01200 

0073 

GO  TO  202 

SFN012I0 

0074 

110  CONTINUE 

SEN01270 

0075 

DO  120  1‘ltNPAR 

SENOI  2  30 

0076 

X2H( 1 1  •  PARI It 

SEN01240 

0077 

120  CONTINUE 

SEN01250 

0078 

00  200  J-liNPAR 

SFN01260 

0074 

IFINEXOPS.NE.OI  GO  TO  IIS 

SEN01270 

0080 

1FIKTEST1JI.EQ.0I  GO  TO  200 

SENOI 280 

0081 

115  IFINEXOP4.EQ.OI  GO  TO  121 

SFN01290 

0082 

CALL  LNULT 1 OtOELNUt  DEM. OUI 

SENOIIOO 

C  COMPUTE  OIOEL  PI/OAIJI  ANO  DIFI/OAIJI 

USING  CENTRAL 

OIFFERENCING. 

SENOI 310 

008  3 

121  PARI J 1  •PARIJI  »  DPARIJI 

SEN01320 

0084 

CALL  RESTNT 40. 0F1 

SEN01330 

0085 

1FIM.EQ.0I  GO  TO  126 

SEN01340 

0086 

00  125  1»1,M 

SEN01350 

0087 

CALL  RESTNTI ItRJl || 1 

SENOI 360 

0088 

ifirjiii.gt.avali  go  to  125 

SEN0I370 

0084 

123  OPARUI-0.1*OPAR(J| 

SENOI 380 

0040 

PARIJI  »  X2HIJI 

SENOI 390 

0091 

MRITEI6.124I  J, DPARIJI 

SEN01400 

0092 

124  FORMAT  1 16H  RESETTING  DPARt»I2.3H|6 

.G14.5I 

SEN01410 

0043 

IFIDPARI JI.EO.O. 1  CO  TO  201 

SEN01420 

0094 

GO  TO  121 

SEN01430 

0095 

125  CONTINUE 

SEN01440 

0096 

126  IFIM2.EQ.0I  GO  TO  128 

SEN01450 

0097 

00  127  l>l,HZ 

SEN01460 

0098 

M2 PI >9*1 

SEN01470 

0099 

CALL  RESTNTINZPItRJINZPIII 

SENOI 480 

0100 

127  CONTINUE 

SEN01490 

0101 

128  IFINEX0P4.EQ.0I  GO  TO  129 

SENOI 500 

0102 

CALL  LNULTIltOELMU.DEMtDUl 

SEN01510 

0103 

129  CALL  GRADI2I 

SEN01520 

0104 

DO  130  l-l.N 

SENOI 530 

0105 

DELXIII-DELXOII  1 

SEN01540 

0106 

130  CONTINUE 

SEN01550 

0107 

0EM>2.0*DPARIJI 

SENOI 560 

0108 

PAR  1 J laPARI J 1  -  DEM 

SEN015T0 

0109 

CALL  RESTNT10. VAL 1 

SEN015B0 

0110 

IFIM.Ed.0l  GO  TO  136 

SEN01590 

0111 

DO  135  l-ltN 

SEN01600 

0112 

CALL  RESTNT 1 ItRJl III 

SENOI 6 10 

0113 

IFIRJIII.GT.AVALI  GO  TO  135 

SEN01620 

0114 

GO  tO  123 

SEN01630 

0115 

135  CONTINUE 

SEN01640 

0116 

136  IFIM2.EQ.0I  GO  TO  138 

SEN01650 

0117 

DO  137  l>ltNZ 

SEN01660 

0118 

N2PI«M*I 

SEN01670 

0119 

CALL  RESTNT IM2P I . RJI N2P1 1 1 

SEN01680 

0120 

137  CONTINUE 

SEN01690 

0121 

138  IFINEXOP4.EQ.OI  GO  TO  139 

SENOI700 

0122 

CALL  LMULTI2.0ELMU.0EM.0UI 

SEN01710 

0123 

139  CALL  GRA0I2I 

SENOI  720 

0124 

DF-IDF-VALI/DEM 

SENOI730 

0125 

DO  140  l«ltN 

SENOI 740 

0126 

OELX 1 1 1*1  DELXI I 1  -  DELXOI 1 1  I/OEM 

SEN01T50 

0127 

140  CONTINUE 

SENOI TOO 

0128 

PARIJI  ■  X2HIJI 

SENOI TTO 

C  HAVING  ALREADY  FACTORED  A.  SOLVE  A6X-8  FOR  X.  SENOI 7RO 

C  WHERE  6-OELX  AND  X«DX/OAIJ|.  SENOI T90 
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V 

V 

►  ' 


.  1 

0129 

CALL  1NVEKSI2I 

1 

C  PRINT  007  01/08(03 

0130 

NR1TEI6.150I  J 

0131 

150 

FOHNATT47M  X-0ER1VAT 1 VES  ARE  KITH  RESPECT 

TO  PARAMETER 

r 

0132 

DO  170  I ■ 1 ,N»  6 

0133 

U-M1N0I  I«S,NI 

0134 

HRITEI6.160I  I JJ.OELXTJJ) ,  JJ-I.III 

0133 

160 

FORMAT (  614H  OXt , 12,2HI-,Gi4.71  1 

0136 

170 

CONTINUE 

C  CALCULATION  OF  CP  VARIABLE  SENSITtVITV 

0137 

IFINEX0P6.NE.il  CO  TO  670 

0138 

I  NOE  X-l 

0139 

CALL  TRANSI1N0EXI 

0140 

NR1TEI6.610)  J 

0141 

410 

FQRMATI2X, •CORRESPONDING  DERIVATIVES  OF  GP 

VARIABLES  HAT 

U*  PARANE  TER  * »  I  2i  ‘ARE*  1 

0142 

DO  680  I>l.N>6 

0143 

I1»H1N0TI*5.N) 

| 

0144 

KRITEI6, 16011 JJ.OELXTJJ). JJ- 1,113 

0145 

680 

CONTINUE 

1 

0146 

00  690  1*1  »N 

1 

0147 

690 

D£LX1I)«TDELX(II 

0148 

670 

CONTINUE 

0149 

IFINEX0P4.E9.0I  GO  TO  375 

0150 

IF (N.EO.OI  GO  TO  301 

0151 

CALL  LNULT 1 3. DEL MU. DEN. DU 1 

0152 

MR1TE46. 3501  J 

0153 

350 

FORMAT  1/  41H  U-0ERIVAT1VES  KITH  RESPECT  TO  PARAMETER 

0154 

00  351  I'liK.A 

0155 

11  -  MINOT  1  *5,  Ml 

01)6 

KRITE 16.3521  1 1 JJ.DUt J J 1 1 ,JJ«I , II I 

01)7 

352 

FORMAT I6I4H  DUI ■ 1 2. 2HI« >614. 71 1 

i 

0158 

351 

CONTINUE 

01)9 

301 

1FIM2.E9.0I  GO  TO  375 

0160 

CALL  LNULTI4,DELNU.0EM.DUI 

0161 

HR1TET6. 3601  J 

0162 

360 

FORMAT  1 /  41H  W-0ERIVAT1  VES  KITH  RESPECT  TO  PARAMETER 

0163 

DO  361  ■•ltMEtA 

0164 

II  «  MINOT  1*5,621 

i 

0165 

MRITEI4.362I  T I  JJ.DUT JJ*NI),  JJ-I.IO 

0166 

362 

FORMAT I6I4H  DMT , I2.2MI- ,614.71 1 

i 

0147 

341 

CONTINUE 

0168 

375 

CONTINUE 

C  COMPUTE  OF/ DAI J 1 • 

| 

0169 

00  ISO  I»  1 »  N 

I 

0170 

OF  •  OF  *  X3MT  |  l*DELXT  1 1 

0171 

ISO 

CONTINUE 

C  PRINT  Of /DAI  31  • 

0172 

KR1TEI 6, 1901  DF 

l 

0173 

190 

FORMAT 1/  10X,13HDFIX|K)1/DA>  ,C 14.4/10H  **< 

t 

0174 

200 

CONTINUE 

0175 

GO  TO  202 

i 

0176 

201 

KRITE 1 6. 106)  J 

i 

0177 

106 

FORMAT  T IM  . 2 1HTE ANIMATING  PARAMETER, 1 3, 16M 

DUE  TO  OPAR  - 

j 

0178 

GO  TO  200 

0179 

202 

NT 3*NP 

0180 

CALL  REJECT 

0181 

DO  205  l>l,N 

0182 

DEL IT  1 1  *  OELTXTII 

0183 

205 

DELXOTII  •  DELTTII 

0184 

400 

CONTINUE 

, 

C185 

500 

RETURN 

0184 

f  NO 

Figure  34. — continued 
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SENOISOO 

SEN01810 

SEN01820 

SENOIAIO 

SEN0IS40 

scxoiiw 

SEN01860 

send is  to 
senoiaao 

SEND  1490 
SEN0I900 
SEN01910 
SEN01920 
SEN0I930 
SEN01943 
SENOI950 
SF NO I960 
SEN0I970 
SENOI980 
SEN01990 
SEM02000 
SEN02010 
SE *02020 
SEN02030 
SEN02040 
SEN02050 
SEN02040 
SEN020T0 
SEN02080 
SEN02090 
SEN02100 
SEN02110 
SEN02120 
SEN02130 
SEN02140 
SEN02150 
SEN02I60 
SEN021T0 
SEN021S0 
SEN02190 
SEN02200 
SEN02210 
SEN02220 
SEN02230 
SEN02240 
SEN02250 
SEN02260 
SEN022T0 
SEN022SG 
SEN02290 
SEN02300 
SEN02310 
SEN02320 
SEN02330 
SEN0234G 
SEN023S0 
SE  NO 2 340 
SEN02320 
SEN02300 
SEN0239Q 
SEN02400 
SEN02410 
SEN02420 
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FORTRAN  IV  6  LEVEL  21 


S£T 


OAT 6  •  802*2 


12/14/52 


0001 

c 

SUBROUTINE  SET  ITNNAX1 

c 

£ 

FEOUART  1971 

c 

SET  STORES  THE  TINE  AT  WHICH  THE 

<7002 

CONNON  /TSW /  NSMH 

0003 

r 

CONNON  /TNX/  TNQ.EXT .EKT90 

V 

C 

£ 

SECOND  GIVES  JOS  CFU  EXECUTION 

0004 

CALL  SECOND  ITNOI 

0003 

EXI*TMNAX»TNO 

000* 

EXT90"  TNO  *  0.90*TPNAX 

0007 

NSWH-1 

0000 

RETURN 

0000 

ENO 

SET 00040 
SET00050 
SET00060 
SET00070 
SETOOOSO 
SET 00090 
SETOOIOO 
SE T001 10 
SET00I20 
SET00130 
SET00140 
SET00130 
SETOOIfrO 
SET00170 
SET001A0 
SETOOIOO 


Figure  35. — Subroutine  SET. 


6.28  ST0RE 

Subroutine  ST0RE  stores  the  values  of  the  current  point  x  and 
the  associated  values  of  F,  w,  and  G  in  a  temporary  storage  area.  These 
values  are  restored  by  a  call  to  REJECT.  The  listing  for  subroutine 
ST0RE  appears  as  Figure  36. 


FOkfKAN 

IV  G  LEVEL  21  STORE  DATE  ■ 

S02A2 

12/13721 

0001 

SUBROUTINE  STORE 

STOOOOAO 

C 

ST000050 

c 

OCTOBER  1970 

ST00006Q 

c 

ST000070 

c 

STORE  STORES  THE  VALUES  OF  THE  CURRENT  POINT  AND 

THE  ASSOCIATED 

STOOOOSO 

c 

VALUES  OF  THE  FUNCTIONS  IN  A  TERROR ARV  AREA. 

ST000090 

0002 

IMPLICIT  REAL*SIA-H.O-ZI 

ST000100 

000) 

REAL *4  RHOIN. RATIO. ERSI.THETAO 

ST0001 10 

0004 

CONMON/SHAAE/Xl 45 1 ,OELI  AS  1 ,A( 45,431 ,N,N, NN.NF 

l.NNI 

ST000120 

0003 

CONNON  /EOAL/  H.  HI.  N2 

ST000130 

0006 

CUNHUN/VALUE/F, G, PO, RSIGMA.RJ 1 90) .RHO 

ST000140 

000 1 

CONNON/CAST/OEL  XI 451 .DEL X 01 451. RH01 N. RATIO, ERS1 .THETAO, 

ST000150 

iRSIGl,Gl.Xll45I.X2(45l,X3l45l ,XR2 ( 43) ,XRl 1451 

.PR1. 

STOOOIAO 

2RR2.R1.F1,RJ1I90),D0TT.PGRAD|ASI,0IAGIASI. 

ST0001 TO 

3  PREV3, AOELX ,  NICER.  NUNINI.  NPHASE.  NSATIS 

stoooibo 

0008 

DO  10  l«l,N 

ST00U.90 

00  OV 

10 

xmi*xii> 

ST000200 

ooio 

NN2>N*N2 

ST000210 

ooii 

00  20  >l,NN2 

ST000220 

0012 

20 

RJIUI-RJU) 

ST000230 

0013 

Pl-PO 

ST000240 

0014 

Fl-F 

ST000230 

0013 

Gl-G 

ST 0002*0 

001* 

RSICfRSIGMA 

ST000270 

0017 

HI  aN 

ST0002RO 

0018 

RETURN 

STOOD 290 

0019 

ENO 

ST000300 

Figure  36. — Subroutine  ST0RE. 
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6.29  TCHECK 

If  the  print  option  has  been  set  so  that  printouts  occur  only 
after  a  subproblem  has  been  called,  then  TCHECK  is  called  after  each 
iteration  of  the  algorithm  used  to  minimize  the  W  function.  The  elapsed 
time  is  computed  as  is  done  in  TIMEC  but  it  is  not  printed  out.  If  the 
elapsed  time  exceeds  90  percent  of  the  estimated  time  limit,  then  the 
print  option  is  changed  so  a  printout  occurs  after  each  iteration  of  the 
procedure  being  used  to  minimize  the  W  function.  If  the  elapsed  time 
exceeds  the  user  specified  time  limit,  then  the  routine  will  cause  ter¬ 
mination  of  the  attempt  to  solve  the  problem  by  setting  NSWW  equal  to  2. 
Subroutine  TCHECK  appears  as  Figure  37. 


FORTRAN  IV  G  LEVEL  21 


TCHECK  DATE  •  80243  10/40/07 


0001 


0002 
000  J 

0004 

0004 

0006 


0007 

OOOll 

0004 

0010 

oou 

0012 

0011 

0014 


SUHHOUTINC  TCHECK  TCH00940 

C  FEBUARV  1971  TCH00060 

C  TCH00070 

C  TCHECK  CHECKS  THE  NUMBER  OF  SECONDS  THAT  HAVE  ELAPSED  SINCE  THE  START  ICHOOOR'J 

C  UE  THE  PROBLEM.  IF  THE  SOLUTION  IS  TAKING  LONGER  THAN  90  PER-CENT  TCH00090 

C  OF  THE  ESTIMATED  MAXIMUM  TIME,  A  SNITCH  IS  SET  TO  GIVE  MORE  OUTPUT.  TCH00100 

COMMON  /TSM/  NSWM  TCHOOUO 

COMMON  /OPTNS/  NTI,NT2,NT3,NT4,NT5,NT6,NT7,NTS,NI9,NT10  TCH00I2O 

COMMON  /TMX/  TMOiEAT ,EXT90  TCHQ0139 

CALL  SECOND  (SECSI  TCH0014O 

IF  I SECS.LT.EXT90I  RETURN  TCH00140 

C  GETTING  CLOSE  TO  EXCEEDING  THE  TINE  LIMIT  SET  OUTPUT  OPTION  TO  GIVE  TCH0016J 

C  MORE  OUTPUT.  TCHO0I70 

NT3»2  TCHOOlftD 

X*  SECS  -  TMO  TCHOOl'-O 

NR1 7EI6, 101  X  TCN00200 

10  FORMAT  (ftX,4HTIME«.F9.3,8H  SECONDS  I  TCHJ0210 

IF  |  SECS  .GT.  EXT  I  NSNM-Z  TCH03220 

CALL  OUTPUT  111  TCH0U233 

RETURN  TCH00240 

ENO  TCH00240 


Figure  37. — Subroutine  TCHECK. 


6.30  TIMEC 


Subroutine  TIMEC  calls  the  same  routine  as  called  by  SET  to  ob¬ 
tain  the  current  status  of  the  computer's  time  clock.  It  then  prints 
out  the  elapsed  time  since  the  call  to  SET.  If  the  elapsed  time  is 
greater  than  the  user  specified  maximum  time  limit,  TIMEC  will  cause 
the  termination  of  the  attempt  to  solve  the  problem  by  setting  NSWW 
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equal  to  2.  The  listing  for  subroutine  TIMEC  is  shown  below  in  Figure 


38. 


FORTRAN  IV  C  LEVEL  21  TIMEC  DATE  ■  80242  12/15/S} 


0001 

SUBROUTINE  TIMEC 

TIMOOOSO 

c 

TIM00040 

c 

FEBUARY  1971 

TIMOOOSO 

c 

TIN00060 

c 

TIMEC  CHECKS  THE  NUMBER  OF  SECONDS  THAT  HAVE  ELAPSED  SINCE 

THE  START 

iiMoooro 

t 

OF  THE  PROBLEM.  IT  PRINTS  THIS  NUMBER.  IF  THE  SOLUTION  IS  i 

TAKING 

TIMOOOSO 

c 

LUNGER  THAN  THE  ESTIMATED  MAXIMUM  TIME*  A  SWITCH  IS 

SET  TO 

TERMINATE 

TIMU0090 

c 

THE  NUN. 

TIM00100 

0002 

COMMON  /TSM/  NSUH 

TIM001I0 

0002 

COMMON  /TNX/  TNO,EXT,EXT90 

TIM00120 

c 

TINOOISO 

c 

SECOND  GIVES  008  CPJ  EXECUTION  TIME  IN  1/1000  O'  A 

SECOND 

UH00I40 

c 

TIN001S0 

0004 

CALL  SECOND  ISECSI 

T IN00160 

000% 

X»SECS-TN0 

TIH00170 

0006 

WRITE  16,201  X 

T1M00180 

000/ 

IF  1  SECS. LT. EXT  1  GO  TO  10 

TIM00190 

oooa 

NSaM>2 

TIM00200 

ooov 

10 

RETURN 

TIN00210 

c 

TIM00220 

0010 

20 

FORMAT  I6X,SHTINE-F9.1,8M  SECONDS! 

TIM002S0 

0011 

ENO 

TIM00240 

Figure  38. — Subroutine  TIMEC. 

6.31  TRANS 

TRANS  is  a  special  purpose  subroutine  that  performs  an  exponen¬ 
tial  transformation  of  the  solution  vector  and  sensitivity  results. 

This  subroutine  is  useful  when  solving  problem  c(x)  ,  the  convex 

equivalent  of  the  geometric  programming  problem  G(t)  .  The  former 

-x . 

problem  is  obtained  from  the  latter  via  the  transformation  t  =  e  , 
i=l,...,n.  TRANS  is  used  to  transform  the  results  to  the  original  space 
of  the  variable  t  .  The  motivation  for  coding  this  subroutine  was  the 
computational  difficulty  experienced  by  some  users  while  solving  geo¬ 
metric  programming  problems  with  SUMT.  TRANS  is  called  by  the  subrou¬ 
tines  0UTPUT  and  SENS.  Figure  39  is  a  listing  of  subroutine  TRANS. 


6 . 32  XM0VE 


Subroutine  XM0VE  contains  most  of  the  logic  of  the  algorithm  used 


to  minimize  the  W  function  for  a  given  value  of  r 


(RH0) .  First,  a 
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PONT RAN  IV  G  I  IVK  21  TKANS  OAT E  •  802*1  12/19/10 


0001 

SUiUUJTINE  TALNSI INDEX) 

TRA00350 

c 

SUBROUTINE  TKANS  CUJtu  BV  GHAENMI9/9)  TRANSFORMS  THE  OPTIMAL 

TJUQOOtv) 

CA*  AN)  SENSITIVITY  K  t  SUL  T  ST  0  T  SPALL  VMERE 

f  R4000  /0 

L 

I»E/<M-X|.  THIS  BOUTIN*  IS  USFFJLl  WHEN  OSIN  SOLVES  THE  CONVEX 

TR40JOPD 

c 

fDJIVAllCNT  or  GP  PROBLEMS  OBTAINED  or  THE  ABOVE  EXPONENTIAL 

TftAOO  .>9  > 

c 

TRANSFORMATION. 

TRAOOIO  1 

0002 

IMPLICIT  NLAL*BIA-H.O-2) 

TRA00110 

0003 

LONH JN/SHAPL/X 1*3) .DELIAS  1 .AIAS.ASI ,N. N . NN, NP1 ,NN1 

TKA3412J 

000* 

COHNON/CNS T/OLL X 1 *3 1 .OELXOIAS It KH0IN.NATI0.EPS 1 .THETAO. 
IKSlul.Gt.XllASI ,X2I*S).XJ1*5I |XR?I*3I,XKII*S) ,PR1, 
2PK2,P1.F1,RJ1I9U).00T/»PGRADI*5>.D1AGI*5). 

3  PKEV1.A0ELX,  NIC TR .  NON1NI.  NPHASE.  NSATIS 

IRA00130 

FRA00I4') 

FRAO0190 

TRAO0I60 

0003 

COMMOm/AdGS/T  X I  *31 , TOEL  XI  *51 

TRA00I70 

0006 

C0HMUN/ABG6/XXI  *51 

TRAOOIRO 

000/ 

IFIInJtx.EJ.il  GO  TO  20 

TRA001V.') 

000* 

DU  10  1*1. N 

TKA0020O 

0000 

IXII  l*XIII 

TRA00210 

0010 

xii)*jexpi-xid) 

TRA00220 

0011 

to  x>m*xm 

I RAOO? JO 

0012 

RE  TORN 

TRA1J2A0 

001 1 

20  DO  30  1*1, H 

TRA00250 

ooi* 

TOELXI 1 l*DEL  XI 1 1 

FRA J026J 

0013 

10  otLxl 1 |*-DEXP(-XI I ) )*DELXI I ) 

TRA00270 

0016 

RETURN 

T  RAO 0280 

001/ 

mo 

TRA00290 

Figure  39. — Subroutine  TRANS. 


direction  is  chosen  in  which  it  is  believed  the  W  function  will  initial¬ 
ly  decrease.  Then  0PT  is  used  to  find  a  minimum  of  the  w  function  in 
that  direction.  Having  generated  a  new  point  giving  a  lower  value  of 
the  W  function,  XM0VE  returns  control  to  B0DY. 

Subroutine  XM0VE  contains  three  procedures  for  generating  a  di¬ 
rection  vector.  The  choice  of  the  procedure  used  ic  controlled  by  ex¬ 
perimental  option  2  (NEX0P2).  If  NEX0P2  =  1,  then  Newton's  method  is 

used  to  choose  the  direction  vector  S  ,  that  S1  -  -[V2W(x1,rk)]_1VW(x1) 

If  V^C  x^,r^)  is  not  a  positive  definite  matrix,  S*  is  generated  by 
rules  outlined  on  page  167  in  Fiacco  and  McCormick  [10]. 

If  NEX0P2  =  2,  then  the  negative  of  the  gradient  is  used  as  the 
i  i  k 

direction  vector,  S  =  -VW(x  ,r  )  . 

When  NEX0P2  =  3  a  variable  matrix  method  is  used  to  generate 

S1  .  This  method  is  McCormick's  modification  of  the  Fletcher-Power- 
Davidon,  as  described  on  pages  170-175  in  Fiacco  and  McCormick  [10]. 
Subroutine  XM0VE  is  shown  as  Figure  40. 
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FORTRAN  IV  &  LEVEL  21 


XNQVE  DATE  •  802A2  12/16/26 


0001 

SUttkUUTINE  XMOVE 

KM000040 

C 

X MOO 0050 

C 

MARCH  1971 

XM000060 

C 

XN00007J 

C  XMOVE  DETERMINES  THE  VECTOR  ALONG  WHICH  THE  SEARCH 

FOR  A 

MINIMUM  IS 

XN000080 

C  USING  OPT. 

XH030390 

0002 

IMPLICIT  AEAL-8IA-H.0-ZI 

XMOOOIOO 

000) 

REAL94  RHOIN. RATIO, EPSI .THETAO. XEP I , XEP2 

XMOOOUO 

0004 

C0HM0N/SHARE/Xt45),DELt45l,At45.45t,N,M,NN,NPl 

,  NMl 

XMOOO 120 

0005 

CUMMUN/CAST/0ELX(45t iDELXOI 45) .RHOIN, RAT  10. EPSI. THETAO. 

XN000130 

1RS1G1.G1.X1145I.X2I45I.X3I45I .XR2I45) .XR1I45I.I 

PRl, 

XMOOO 140 

2PR2.P1.F1.RJ1I90),D0TT.PGRADI45).DIAG:4SI. 

XMO00I50 

3  PREV3.A0ELX,  NTCTR,  NUMINI.  NPHASE.  NSATIS 

XM000160 

0006 

COMMON/E  XPOPT/ NEXOPl .NEX0P2. NE  XOP). NE  X0P4.NEX0P  5. XEP 1 

I.XEP2 

XM0001 70 

0007 

COMMON/XVE/ SIGi 451. V VI 45) .XXXI 45) .OCLL 1451 

KM000180 

C-- 

NEX0P2  DETERMINES  HOW  MOVE  IS  TO  BE  MADE 

XMOOOIOO 

C  NEX0P2  •  1  USE  NODIMEO  NEWTON  RAPHSON  NETHOO. 

XM00020L 

C 

«  2  USE  MOOIF1EO  NEKTON  RAPHSON  METHOD.  BUT  ROD 

DELXO  TO 

XM0002IO 

c 

ORTHOGONAL  MOVE  VECTOR  IF  HESSIAN  IS  INDEFINITE. 

XM000220 

c 

•  3  USE  STEEPEST  DESCENT  METHOD. 

XM000230 

c 

•  4  USE  MCCORMICK. S  MODIFICATION  OF  THE  FLETCHER- 

POWELL 

XM000240 

c 

NETHOO. 

XM000250 

000 8 

GO  TO  110.10,130,30),  NEX0P2 

XM000260 

c— 

NEKTON  -RAPH  KITH  KHAT EVER  NETHOO  IS  IN  INVERSE 

XM000270 

0009 

10 

CALL  GRAO  111 

XM000280 

c— 

ONE  111  MEANS  ACCUMULATE  MATRIX  OF  SECOND  PARTIAL 

DERIVATIVES 

XM000290 

0010 

CALL  SECORD  11) 

XM000300 

0011 

00  20  I-l.N 

XMOOO) 10 

0012 

20 

OELXII I.OELXOI  1 1 

XMOOO370 

0013 

CALL  INVERS  111 

XM000330 

C  IF  A  NONPOSITVE  PIVOT  IS  ENCOUNTERED  IN  INVERSE  AN 

ATTEMPT  IS  HADE  TO 

X MO 00340 

C  COMPUTE  A  VECTOR  HAVING  A  POSITIVE  DOT  PRODUCT  MITH  A  NEGATIVE 

XMOOO 350 

C  EIGENVECTOR  ANO  THE  NEGATIVE  OF  DEL  P. 

XMOOO 360 

0014 

CALL  STORE 

XM030370 

0015 

CALL  OPT 

XN000380 

0016 

RETURN 

XM000390 

C— 

F-P-O-NCC  MOVE 

XM000400 

0017 

30 

CALL  GRAO  121 

XN000410 

C — 

NN  IS  NO.  OF  MOVES  FOR  THIS  VALUE  OF  RHO 

XM000420 

0018 

IF  INN.NE.OI  GO  TO  70 

XN000430 

0019 

40 

IREP«0 

KM000440 

0020 

IT-0 

XM000450 

C  — 

SET  INITIAL  GUESS  INVERS  MATRIX  OF  SECOND  PARTIAL 

DERIVATIVES 

XM000460 

c— 

use  partial  inverse  if  known 

XM000470 

0021 

00  50  1-1, N 

XM000480 

0022 

00  50  J-l.N 

XN000490 

0023 

50 

Al  1. JI-0.0 

XM000500 

0024 

00  60  1-1 .N 

XN000510 

0025 

60 

Al 1,11-1.0 

XM000520 

0026 

TO 

DO  80  I-l.N 

XM000530 

0027 

80 

DELXI 1 l-OELXOI 1 1 

XM000540 

0028 

IF  IIREP.GT.N)  GO  TO  40 

XM000550 

0029 

IF  IIT.EO.OI  GO  TO  130 

XM000560 

0030 

DO  90  1-1, N 

XM000570 

0031 

sight -x  ii>- xx  xiii 

XN000580 

0032 

90 

VVII)-DELLI2 l-OELXOI ! I 

XMOOO 590 

C— 

NEGATIVE  GRADIENT  STORED  ANO  CONFUTED 

XM000600 

C— 

COMPUTE  NT 

XN000610 

003) 

DO  100  1-ltN 

XM000620 

00)4 

OELXII 1-0.0 

XM000630 

0035 

DO  100  3-1, N 

XM000640 

0036 

100 

OELXI ll-OELXIl l-AII.JI-YYI  J> 

XM000650 

c— 

COMPUTE  VISIG  —MV  1-1 

XM000660 

00)7 

2CON-0.0 

XM000670 

00)8 

DO  110  l-l.N 

XM000680 

00)9 

110 

2C0N-ZC0N4rrill«ISiGII  l-DELXl  III 

XN000690 

0040 

IF  I2CON.EO.O.O)  GO  TO  1)0 

XMOOO 700 

0041 

IREP-IREP-1 

XN0007IO 

0042 

2C-1./2CON 

XM000720 

c— 

UPDATE  H  MATRIX  USING  NCC  FORMULA  MIEN  SCALAR 

NE  0 

XM000730 

004) 

00  120  I-l.N 

XMOOO  740 

0044 

T1-2C*ISIGII l-DELXIII) 

XN000750 

0045 

00  120  2-I.N 

XM000760 

0046 

Al  1,  JI-AII,  JI»T1»I-DELXIJ)*SIGUII 

XMOOO 7  70 

0047 

120 

Al  J.II-AI I.JI 

XN000T80 

C— 

STORE  CURRENT  POINT  ANO  CURRENT  GRAOIENT  INEGI 

XMOOO 790 

0048 

1)0 

00  140  I-l.N 

XM000800 

0049 

XXXID-XIII 

XMOOOSIO 

0090 

140 

OELLII l-OELXOI II 

XM000820 

0051 

00  150  I-l.N 

XN000830 

0092 

OELXII 1-0.0 

XM000840 

009) 

00  190  J-l.N 

XMOOO 850 

0094 

150 

OELXII l-OELXI 1 1 -All , Jt *DELXO t  Jl 

XN000860 

0095 

2C1-0.0 

XMOOO* 70 

0096 

00  160  I-l.N 

XM000880 

0097 

160 

2C  l-OELXI 1 1-62-2C1 

XM000890 

0046 

7CI-OSORTI2l.il 

XM000900 

Figure  40. — Subroutine  XM0VE 
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00*9 

00  170  1*1, N 

X *000910 

0000 

170 

0£LXII|*DELXIII/2C1 

**000920 

0061 

CALL  STORE 

X MO 009 50 

0062 

CALL  OPT 

XMQ00940 

006* 

1 T- *  T*1 

XH0009S0 

0064 

RETURN 

XM000960 

006* 

180 

CONTINUE 

X* 000970 

C  STEEPEST  OESCENT 

XM0009S0 

0066 

CALL  GRAD  121 

X *000990 

006  7 

00  190  1*1, N 

X *001000 

0068 

190 

OELXtll*DELXO<ll 

X*00!010 

0069 

CALL  STORE 

X*00(0?0 

0070 

CALL  OPT 

X *001 030 

0071 

RETURN 

X *001 090 

0072 

CNO 

X *00 1050 

Figure  40. — continued 
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